# Rigid formic acid simulation variable T equal 100 variable dtequil equal 3.0 variable dt equal 1.0 # Comments from manual, see pag 35 # 1. Initialization: units real boundary p p p dimension 3 atom_style full neighbor 2.0 bin neigh_modify every 1 delay 0 check yes page 300000 one 5000 # 2. Atom definition: read_data 0nacl0h2o.lammps05 mass 1 12.011 #c mass 2 15.999400 # O mass 3 15.999400 # O mass 4 1.008 mass 5 1.008 timestep ${dtequil} # 3. Force fields: pair_style lj/cut/coul/long 14.0 # I J eps(kcal/mol) sigma(A) #I,J=atom_types pair_coeff 1 1 0.0899 3.727 #c pair_coeff 1 2 0.0918 3.4535 pair_coeff 1 3 0.1615 3.2005 pair_coeff 1 4 0.0208 2.2635 pair_coeff 1 5 0.0464 2.3605 pair_coeff 2 2 0.0937 3.18 #o(c-o) pair_coeff 2 3 0.1649 2.927 pair_coeff 2 4 0.0212 1.99 pair_coeff 2 5 0.0473 2.087 pair_coeff 3 3 0.2902 2.674 #o(c=o) pair_coeff 3 4 0.0373 1.737 pair_coeff 3 5 0.0833 1.834 pair_coeff 4 4 0.0048 0.8 #h(c-h) pair_coeff 4 5 0.0107 0.897 pair_coeff 5 5 0.0239 0.994 #h(oh) bond_style harmonic # N k/2(kcal/mol/A^2) r0(A) bond_coeff 1 450 1.35 bond_coeff 2 450 1.213 bond_coeff 3 450 1.096 bond_coeff 4 450 0.98 angle_style harmonic # N k/2(kcal/mol/rad^2) theta0(degree) K_UB(en/dis^2) rub(A) angle_coeff 1 55 125.5 angle_coeff 2 55 125.1 angle_coeff 3 55 109.4 angle_coeff 4 55 106.1 kspace_style pppm 1.0e-4 velocity all create $T 87287 loop all dist gaussian fix a2 all rigid/npt molecule temp $T $T 100 iso 1.0 1.0 1000 fix_modify a2 energy yes thermo_style custom step temp etotal pe epair ebond eangle press vol f_a2 thermo 1 dump 1 all custom 2000 dump.meltq401 id type x y z mass q compute mytemp all temp compute 4 all pressure mytemp fix cc1 all ave/time 500 1 500 c_mytemp c_4 file tmp.T100 mode scalar run 10000 timestep ${dt} unfix a2 fix a2 all rigid/npt molecule temp $T $T 100 iso 1.0 1.0 1000 run 1000000