variable T equal 298 #variable Tl equal 190 variable p equal 1 variable step equal 5 # timestep (fs) variable runtime equal 4800000 variable dumptrj equal 4800 # frequency in timesteps to dump the trajectory variable out equal 1000 # frequency of writing thermo quantieis to output file variable Tdamp equal ${step}*100 # Thermostat parameter variable pdamp equal ${step}*1000 # Barostat parameter variable f_restart equal ${runtime}/10 units real atom_style full boundary p p p special_bonds charmm bond_style harmonic angle_style harmonic dihedral_style charmm improper_style harmonic pair_style hybrid lj/cut 12.0 sw sw0 yukawa 1.0 10.0 pair_modify tail yes mix arithmetic read_data CG_checked.data replicate 4 4 4 #pair_modify mix arithmetic tail yes # Na=TMA, Me=Methyl, Ar=Aromatic, Et=ethyl, O=ether oxygen # sw file sets Ar and Me with mini methane parameters (Bin). # They are overwritten below with OPLS parameters (Yuqing) #pair_coeff * * sw0 ppotma.sw Me Ar Me Me O Cl mW Me Na pair_coeff * * sw0 ppotma.sw mW Ar Me Na O Cl NULL Ar Ar #pair_coeff * * sw0 ppotma.sw NULL NULL NULL NULL NULL Cl mW NULL Na #pair_coeff * * sw ppotma.sw NULL NULL NULL NULL NULL Cl mW NULL NULL pair_coeff * * sw ppotma.sw mW NULL NULL NULL NULL Cl NULL NULL NULL #Use OPLS parameter for all hydrophobic groups (Ar, Me, and O) pair_coeff 2 2 lj/cut 0.11 3.75 # from opls-ua for benzene pair_coeff 8 8 lj/cut 0.11 3.75 # from opls-ua for benzene pair_coeff 9 9 lj/cut 0.11 3.75 # from opls-ua for benzene #pair_coeff 3 3 lj/cut 0.207 3.775 # from opls-ua for CH3 (C1) pair_coeff 3 3 lj/cut 0.175 3.905 # from opls-ua for CH3 (C2) pair_coeff 5 5 lj/cut 0.17 3.00001 # O-O from amber #pair_coeff 2 3 lj/cut 0.1509 3.7625 # arithmetic Ar-Me for CH3 (C1) pair_coeff 2 3 lj/cut 0.1387 3.8275 # arithmetic Ar-Me for CH3 (C2) pair_coeff 3 8 lj/cut 0.1387 3.8275 # arithmetic Ar-Me for CH3 (C2) pair_coeff 3 9 lj/cut 0.1387 3.8275 # arithmetic Ar-Me for CH3 (C2) pair_coeff 2 5 lj/cut 0.1367 3.375 # arithmetic Ar-O pair_coeff 5 8 lj/cut 0.1367 3.375 # arithmetic Ar-O pair_coeff 5 9 lj/cut 0.1367 3.375 # arithmetic Ar-O #pair_coeff 3 5 lj/cut 0.1876 3.3875 # arithmetic Me-O for CH3 (C1) pair_coeff 3 5 lj/cut 0.1725 3.45251 # arithmetic Me-O for CH3 (C2) pair_coeff 1 2 lj/cut 0.17 3.336 # mW-Ar pair_coeff 1 8 lj/cut 0.17 3.336 # mW-Ar pair_coeff 1 9 lj/cut 0.17 3.336 # mW-Ar #pair_coeff 1 3 lj/cut 0.20 3.536 # mW-Me pair_coeff 1 3 lj/cut 0.215 3.4 # mW-Me #pair_coeff 1 5 lj/cut 0.20 3.536 # mW-O pair_coeff 1 5 lj/cut 0.25 3.036 # mW-O pair_coeff 2 6 lj/cut 0.47 3.61 # Ar-Cl xxxxx pair_coeff 6 8 lj/cut 0.47 3.61 # Ar-Cl pair_coeff 6 9 lj/cut 0.47 3.61 # Ar-Cl pair_coeff 2 4 lj/cut 0.15 4.3 # Ar-Na xxxxx pair_coeff 4 8 lj/cut 0.15 4.3 # Ar-Na pair_coeff 4 9 lj/cut 0.15 4.3 # Ar-Na #pair_coeff 3 6 lj/cut 0.20 3.536 # Me-Cl xxxxx pair_coeff 3 6 lj/cut 0.50 3.61 # Me-Cl xxxxx pair_coeff 3 4 lj/cut 0.20 4.3 # Me-Na xxxxx pair_coeff 5 6 lj/cut 0.043 4.699 #3.9589 # O-Cl xx pair_coeff 4 5 lj/cut 0.25 4.3 # Na-O xx # zero all the pair interactions between atom and type 7 pair_coeff 1 7 none pair_coeff 2 7 none pair_coeff 3 7 none pair_coeff 4 7 none pair_coeff 5 7 none pair_coeff 6 7 none pair_coeff 7 7 none pair_coeff 7 8 none pair_coeff 7 9 none #pair_modify tail yes mix geometric pair_coeff 6 6 yukawa 1107.0 pair_coeff 4 4 yukawa 1107.0 #TMA-Ar bond #bond_coeff 10 50.0 2.55 # TMA-Ar #dihedral angle that spans ether oxygen #dihedral_coeff 21 1.80 3 180 0.0 # X CA OS X #control angle offset of TMA to plane of ring #improper_coeff 10 10.0 180 #TMA-CA-CA-CA #group hyd type 1 4 8 10 #delete_bonds hyd multi remove any #delete_atoms group hyd #Coarse-Grain masses #mass 3 15.0347 #mass 7 18.015 #mass 9 74.144 #the rings and attached methyl groups are rigid. #for the full polymer you will have to write a script to set each separate ring as different molecules group water type 1 group tman type 4 group tma type 4 group cl type 6 group ring1 type 8 9 group ring2 type 2 group notring subtract all ring1 ring2 compute 1 water msd compute 2 tman msd compute 3 cl msd compute petp tma pe/atom pair compute petmap tma reduce sum c_petp thermo_style custom step temp epair ebond eangle pe etotal press vol lx ly lz c_petmap c_1[4] c_2[4] c_3[4] edihed eimp emol thermo_modify flush yes minimize 0.01 0.001 100 100 minimize 0.01 0.001 100 100 #fix 1 all shake .0001 20 0 b 4 6 a 4 velocity all create $T 1402516963 mom yes rot yes dist gaussian loop local #Make sure NVE conserves energy in the full polymer for this timestep timestep ${step} #use rigid/npt for rings. Everything else is normal NPT fix 1 all momentum 1 linear 1 1 1 #fix 4 all npt temp $T $T ${Tdamp} iso $p $p ${pdamp} fix 2 notring npt temp $T $T ${Tdamp} iso $p $p ${pdamp} dilate all fix 3 ring1 rigid/nvt/small molecule temp $T $T ${Tdamp} fix 4 ring2 rigid/nvt/small molecule temp $T $T ${Tdamp} dump 1 all custom ${dumptrj} traj.lammpstrj id type xs ys zs #ix iy iz thermo ${out} run ${runtime} write_data out.data