#!usr/bin/python
import sys
from random import choice, random
import os
import subprocess

N = 10
ndenmax = 10

   
if type(N) != int:
    print 'N must be an integer.'
   
Ntag = str(N)
ntag = str(ndenmax)
xtag = str(xmax)
ytag = str(xmax)
ztag = str(xmax)
eptag = str(xmax)
sigtag = str(xmax)
qtag = str(xmax)

top_path = ('materials' + '_' + Ntag + '.' + ntag + '_' + xtag + '.' + ytag
   		 + '.' + ztag + '_' + eptag + '.' + sigtag + '_' + qtag)

if not os.path.exists(top_path):
    os.mkdir(top_path) 
def drange(start, stop, step):
    r = start
    while r < stop:
    yield r
    r+= step

nden0 = drange(1, ndenmax, ndenp)
ndendim = [nden for nden in nden0]

x0 = drange(0, xmax + xp, xp)
xdim = [x for x in x0]
    
y0 = drange(0, ymax + yp, yp)
ydim = [y for y in y0]
    
z0 = drange(0, zmax + zp, zp)
zdim = [z for z in z0]
    
ep0 = drange(0, epmax + epp, epp)
epdim = [ep for ep in ep0]
   
sig0 = drange(0, sigmax + sigp, sigp)
sigdim = [sig for sig in sig0]
    
    #open mat_stats.txt, to track material data    
mat_stats = open(os.path.abspath(top_path)+ '/mat_stats.txt', 'w')
mat_stat_heading = ('\nBOUNDARIES\nNumber of particles: ' + Ntag +
                  	'\nnumber density:   ' + ntag + '\nx-coordinate: ' +
			xtag + '\ny-coordinate: ' + ytag + '\nz-coordinate: ' +
			 ztag + '\nEpsilon: ' + eptag + '\nSigma: ' + sigtag 
			+ '\nCharge: ' + qtag + '\n\n' +
			'#name     number density     xdim     ydim     '+
			'zdim     total particles     net charge\n')
mat_stats.write(mat_stat_heading)
    
#MAT-XXX loop...
for i in range(N + 1):

    mat_name = 'MAT-' + str(i)
        
	#make MAT-XXX directory
    os.mkdir(top_path+'/'+mat_name)
	
	#open .cif file
    cif_file = open(os.path.abspath(top_path) + '/'+mat_name + '/' + 
			mat_name+'.cif', 'w')
	
	#open force_field_mixing_rules.def
    mixing_rules = open(os.path.abspath(top_path) + '/'+mat_name +
			'/force_field_mixing_rules.def', 'w')
        
	#open pseudo_atoms.def
    pseudo_atoms = open(os.path.abspath(top_path) + '/'+mat_name + 
			'/pseudo_atoms.def', 'w')
	
	#open force_field.def
    force_field = open(os.path.abspath(top_path) + '/'+mat_name +
			'/force_field.def', 'w')

    N_ = choice(ndendim)
    n_ = int(N_)
    xdim_ = choice(xdim)
    ydim_ = choice(ydim)
    zdim_ = choice(zdim)
    nden_ = xdim_ * ydim_ * zdim_ / N_
	

#

    cif_heading = ('material' + str(i) + 
			'\n\nloop_\n_symmetry_equiv_pos_as_xyz\n  x,y,z\n'+
			'_cell_length_a          ' + str(xdim_) +
			'\n_cell_length_b          ' + str(ydim_) +
			'\n_cell_length_c          ' + str(zdim_) + 
			'\n_cell_angle_alpha       90.0000\n' +
			'_cell_angle_beta        90.0000\n' +
			'_cell_angle_gamma       90.0000\nloop_\n' +
			'_atom_site_label\n_atom_site_fract_x\n' +
			'_atom_site_fract_y\n_atom_site_fract_z\n'+
			'_atom_site_charge\n')
    cif_file.write(cif_heading)

	
    mixing_heading = ('# general rule for shifted vs truncated\nshifted\n' +
			'# general rule for tailcorrections\nno\n' +
			'# number of defined interactions\n' + str(n_+4) +  #check these + XXX values
			'\n# type interaction\n')
    mixing_rules.write(mixing_heading)
        
	#mixing_output.write(mixing_header)
        
    pseudo_heading = ('#number of pseudo atoms\n' + str(n_+9) + 
			'\n#type          print    as     chem     oxidation' +
			'     mass       charge     polarization     ' +
			'B-factor     radii    connectivity     anisotropic' +
			'   anisotrop-type  tinker-type\n')
    pseudo_atoms.write(pseudo_heading)
       	
    make charges
    q = []
    for k in range(n_):
        q.append(0)
    for l in range(5*(n_)):
        m = choice(range(n_))
        n = choice(range(n_))
        if m == n:
            n = choice(range(n_))
        dq = random() * qmax
        if q[m] + dq <= qmax and q[n] - dq >= -1 * qmax:
            q[m] = float(float(q[m]) + dq)
            q[n] = float(float(q[n]) - dq)
        if q[m] > qmax or q[n] < -1 * qmax:
            q[m] = q[m] - dq
            q[n] = q[n] + dq
    for o in range(5*(n_)):
        m = choice(range(n_))
        n = choice(range(n_))
        if m == n:
            n = choice(range(n_))
        dq = random() * qmax
        if q[m] + dq <= qmax and q[n] - dq >= -1 * qmax:
            q[m] = float(float(q[m]) + dq)
            q[n] = float(float(q[n]) - dq)
        if q[m] > qmax or q[n] < -1 * qmax:
            q[m] = q[m] - dq
            q[n] = q[n] + dq
    p = choice(range(n_))
    q[p] = q[p] - sum(q)
    if sum(q) != 0.000000000000000000000:
        for l in range(5*(n_)):
            m = choice(range(n_))
            n = choice(range(n_))
            if m == n:
                n = choice(range(n_))
            dq = random() * qmax
            if q[m] + dq <= qmax and q[n] - dq >= -1 * qmax:
                q[m] = float(float(q[m]) + dq)
                q[n] = float(float(q[n]) - dq)
            if q[m] > qmax or q[n] < -1 * qmax:
                q[m] = q[m] - dq
                q[n] = q[n] + dq
        for o in range(5*(n_)):
            m = choice(range(n_))
            n = choice(range(n_))
            if m == n:
                n = choice(range(n_))
            dq = random() * qmax
            if q[m] + dq <= qmax and q[n] - dq >= -1 * qmax:
                q[m] = float(float(q[m]) + dq)
                q[n] = float(float(q[n]) - dq)
            if q[m] > qmax or q[n] < -1 * qmax:
                q[m] = q[m] - dq
                q[n] = q[n] + dq
            p = choice(range(n_))
            q[p] = q[p] - sum(q)

    mat_charge = str(sum(q))
    cif_file.write('#NET CHARGE: ' + mat_charge + '\n')
    mat_X_stats = (mat_name + '     ' + str(nden_) + '     ' + str(xdim_) + '     ' + str(ydim_) +
			'     ' + str(zdim_) + '     ' + str(n_) + '     ' + 
			str(sum(q)) + '\n')
    mat_stats.write(mat_X_stats)
	

	#pseudo_atom loop
    for j in range(n_):
        x = choice(range(int(xdim_ + 1)))
        y = choice(range(int(ydim_ + 1)))
        z = choice(range(int(zdim_ + 1)))
        ep = choice(epdim)
        sig = choice(sigdim)
            
        charge = q[j]
        if charge < 0:
           atom_X_cif = ('A' + str(j) + '     ' + str(x) + '     ' +
				str(y) + '     ' + str(z) + '    ' +
				str(charge) + '\n')    
       	cif_file.write(atom_X_cif)
        if charge >= 0:
            atom_X_cif = ('A' + str(j) + '     ' + str(x) + '     ' +
				str(y) + '     ' + str(z) + '     ' +
				str(charge) + '\n')
            cif_file.write(atom_X_cif)
            
        atom_X_mixing = ('A' + str(j) + '          LENNARD_JONES     ' +
				str(sig) + '     ' + str(ep) + '\n')
        mixing_rules.write(atom_X_mixing)
	    
        atom_X_pseudo = ('A' + str(j) + '   yes   C   C   0   12.0   ' +
				str(q[j]) + '   0.0   0.0   1.0  1.00   0   ' +
				'0  absolute   0\n')
        pseudo_atoms.write(atom_X_pseudo)

    adsorbate_mixing = ('CH4_sp3     LENNARD_JONES          148.0  ' +
				'        3.73\nO_co2       LENNARD_JONES    ' +
				'      79.0          3.05\nC_co2       ' +
				'LENNARD_JONES          27.0          2.80\n' +
				'N_n2        LENNARD_JONES          36.0    ' +
				'      3.31\n' +
				'# general mixing rule for Lennard-Jones\n' +
				'Lorentz-Berthelot')
    mixing_rules.write(adsorbate_mixing)

    adsorbate_pseudo = ('CH4_sp3       yes      C      C        0   ' +
				'          16.04246   0.0        0.0        ' +
				'      1.0          1.00     0              ' +
				'  0             absolute        0\nCH3_sp3 ' +
				'      yes      C      C        0           ' +
				'  15.03452   0.0        0.0              ' +
				'1.0          1.00     0                0   ' +
				'          absolute        0\nCH2_sp3       ' +
				'yes      C      C        0             ' +
				'14.02658   0.0        0.0              1.0 ' +
				'         1.00     0                0       ' +
				'      absolute        0\nCH_sp3        yes ' +
				'     C      C        0             13.01864' +
				'   0.0        0.0              1.0         ' +
				' 1.00     0                0             ' +
				'absolute        0\nC_sp3         yes      ' +
				'C      C        0             12.0       ' +
				'0.0        0.0              1.0          ' +
				'1.00     0                0             ' +
				'absolute        0\nC_co2         yes      ' +
				'C      C        0             12.0107    ' +
				'0.7        1.508            1.0          ' +
				'0.720    0                0             ' +
				'absolute        0\nO_co2         yes      ' +
				'O      O        0             15.9994   ' +
				'-0.35       0.9475           1.0          ' +
				'0.68     0                0             ' +
				'relative        0\nN_n2          yes      ' +
				'N      N        0             14.00674  ' +
				'-0.482      0.0              1.0          ' +
				'0.7      0                0             ' +
				'relative        0\nN_com         no       ' +
				'N      N        0             0.0        ' +
				'0.964      0.0              1.0          ' +
				'0.7      0                0             ' +
				'relative        0\n')
    pseudo_atoms.write(adsorbate_pseudo)        
      
    force_field_rules = ('# rules to overwrite\n0\n' +
				'# number of defined interactions\n0\n' +
				'# mixing rules to overwrite\n0')
    force_field.write(force_field_rules)

       
cif_file.close()
mixing_rules.close()
pseudo_atoms.close()
force_field.close()
mat_stats.close()

