# 400 [SDS][Na], 46000 H2O, Al2O3 (9600 Al + 14400 O)
# 164.8566 x 95.18 x 162.50 angstrom
# June 24, 2015
# Paul Huang

# Initialization
units	real
dimension	3
boundary	p	p	p
atom_style	full

read_data	Al2O3-SDS-H2O.data

# In 'real' unit, the time is in the unit of femtosecond, fs.
timestep 	1.0
restart 5000 restart.*


# Structural Definition
#
# Pair Coeffs, 
#
# 1  A - Al of Al2O3
# 2  B - CH2 (C2~C11)
# 3  C - CH2 (C12)
# 4  D - OB
# 5  E - S 
# 6  F - O1~O3
# 7  H - CH3 (C1)
# 8  N - Na
# 9  O - O of Al2O3 
# 10 P - O of H2O
# 11 Q - H of H2O

# Bond Coeffs
#
# 1  A-O
# 2  B-B
# 3  B-C
# 4  B-H
# 5  C-D
# 6  D-E
# 7  E-F
# 8  P-Q

# Angle Coeffs
#
# 1  A-O-A
# 2  B-B-B
# 3  B-B-C
# 4  B-B-H
# 5  B-C-D
# 6  C-D-E
# 7  D-E-F
# 8  F-E-F
# 9  O-A-O
# 10 Q-P-Q

# Dihedral Coeffs
#
# 1  B-B-B-B
# 2  B-B-B-C
# 3  B-B-B-H
# 4  B-B-C-D
# 5  B-C-D-E
# 6  C-D-E-F
# 7  O-A-O-A


# Force field definition 
# LJ Cut0ff = 12.0 angstrom; 
# Coul Couoff = 14.0 angstrom 
# If the distance is larger than 14.0 angstron, 
# the long-range Coul interaction is calculated by ppm method

pair_style	lj/cut/coul/long	12.0	15.0	
kspace_style	pppm	1.0e-5

# pair_coeff epsilon (Kcal/mol), sigma (angstrom) 
# SDS-Na: J. Phys. Chem. B, 2008, 112, 1987-2000.
# Al2O3, ClayFF: J. Phys. Chem. B, 2004, 108, 1255-1266.
# check Al & O !!!
pair_coeff   1	1    0.0000013298  4.7943 
pair_coeff   2	2    0.1180000000  3.9050 
pair_coeff   3	3    0.1180000000  3.9050 
pair_coeff   4	4    0.1700000000  3.0000 
pair_coeff   5	5    0.2500000000  3.5500 
pair_coeff   6	6    0.2000000000  3.1500
pair_coeff   7	7    0.1750000000  3.9050
pair_coeff   8  8    0.1153000000  2.2750
pair_coeff   9  9    0.1554000000  3.5532

pair_coeff   10  10    0.1553625123  3.1660
pair_coeff   11  11    0.0000000000  0.0000


pair_modify mix geometric 


# Bond Coeffs
#
# 1  A-O
# 2  B-B
# 3  B-C
# 4  B-H
# 5  C-D
# 6  D-E
# 7  E-F

bond_style	harmonic

bond_coeff  1    2000.00	1.969
bond_coeff  2     620.00	1.530
bond_coeff  3     620.00	1.530
bond_coeff  4     620.00	1.530
bond_coeff  5     600.00	1.420
bond_coeff  6     600.00	1.580
bond_coeff  7     900.00	1.460
bond_coeff  8     100.00	1.000


# Angle Coeffs
#
# 1  A-O-A
# 2  B-B-B
# 3  B-B-C
# 4  B-B-H
# 5  B-C-D
# 6  C-D-E
# 7  D-E-F
# 8  F-E-F
# 9  O-A-O

angle_style	harmonic

angle_coeff	  1   2000.000	120.29 
angle_coeff	  2    124.300	111.00
angle_coeff	  3    124.300	111.00
angle_coeff	  4    124.300	111.00 
angle_coeff	  5    124.300	109.50 
angle_coeff	  6    124.300	112.60 
angle_coeff	  7    102.000	102.60
angle_coeff	  8    102.000	115.40 
angle_coeff	  9   2000.000	101.12
angle_coeff	 10    300.000	109.47


# Dihedral Coeffs
#
# 1  B-B-B-B
# 2  B-B-B-C
# 3  B-B-B-H
# 4  B-B-C-D
# 5  B-C-D-E
# 6  C-D-E-F
# 7  O-A-O-A

# E = K[1+cos(n*fi-d)]
# Parameter: K, n, d, weighting factor (0,1)

dihedral_style	opls

dihedral_coeff	  1     -4.7135		-3.1350		0.3655		-1.56775		
dihedral_coeff	  2     -4.7135		-3.1350		0.3655		-1.56775
dihedral_coeff	  3     -4.7135		-3.1350		0.3655		-1.56775
dihedral_coeff	  4     -4.7135		-3.1350		0.3655		-1.56775
dihedral_coeff	  5     -4.7135		-3.1350		0.3655		-1.56775
dihedral_coeff	  6     -4.7135		-3.1350		0.3655		-1.56775
dihedral_coeff	  7     -4.7135		-3.1350		0.3655		-1.56775


minimize 1.0e-4 1.0e-6 100 1000
min_style cg

#Define atomic groups

group sds type 2 3 4 5 6 7 
group na type 8
group al2o3 type 1 9
group water type 10 11
group mobile union sds na water

velocity mobile create	300.0	293288	dist gaussian mom yes rot yes

fix 1 al2o3 rigid single force * off off off torque * off off off

# Use the default value
neighbor	0.3	bin
neigh_modify	delay 0	every	1	check	yes	
neigh_modify	exclude molecule al2o3	

# VMD V1.9 with Velocity
#dump  	1   all	custom 100 npt.lammpstrj id type xu yu zu vx vy vz
#dump_modify 1 sort id 

# XYZ
dump  	2   all xyz 1000 structure.*.xyz
dump_modify 2 sort id 

# CFG
# Oxygen of H2O is colored as Nitrogen purposely. 
dump	3 all cfg 1000 cfg.*.cfg mass type xs ys zs vx vy vz fx fy fz
dump_modify 3 element Al C C O S O C Na O N H  
dump_modify 3 sort id

compute tempFree mobile temp
compute pressFree all pressure tempFree

# Output frequcency
thermo 	200
thermo_style	custom	step  vol c_tempFree c_pressFree temp press	
thermo_modify	line	one


# RUN @ NPT 300.0 K, 1atm, 10.0 ns
fix 2 water shake 0.0001 20 0 b 8 a 9
fix	3	mobile	npt	temp	300.0	300.0	100.0	iso	1.0	1.0	1000.0
fix_modify 3 temp tempFree press pressFree

run	5000000