Hello,

I am having an issue with recreating the MARTINI model derived by marrink et. al. in LAMMPS. I have been working on this problem for a while, and have an extensive method developed, but I find that when I turn off the momentum control in the final analysis run, everything freezes and drifts. All units are non-dimensional based by the MARTINI water model. I am trying to simulate a Butane only system (atom type 2 here) with 400 sites and periodic boundary conditions. The non-dimensional mass for butane in this system is set to 0.806538. The same drifting and freezing problem happens when using 400 water sites. Please Help.

Here is the code in my in file:

#Lammps infile water, INNER, W,0, O,400, DPC_,0, MOLECULES,1,OILTYPE

#-----------------------STARTING RUN--------------------------#

units lj

dimension 3

atom_style full

lattice fcc 1.0

read_data atoms.txt

group oil molecule 2

pair_style lj/gromacs/coul/gromacs 1.91489 2.55319 0.0 2.55319

pair_coeff 1 1 1.000000 1.0

pair_coeff 1 2 0.360000 1.0

pair_coeff 1 3 1.000000 1.0

pair_coeff 1 4 1.000000 1.0

pair_coeff 1 5 0.680000 1.0

pair_coeff 1 6 1.000000 1.0

pair_coeff 1 7 1.000000 1.0

pair_coeff 2 2 0.680000 1.0

pair_coeff 2 3 0.360000 1.0

pair_coeff 2 4 0.360000 1.0

pair_coeff 2 5 0.520000 1.0

pair_coeff 2 6 0.360000 1.0

pair_coeff 2 7 0.360000 1.0

pair_coeff 3 3 0.680000 1.0

pair_coeff 3 4 0.680000 1.0

pair_coeff 3 5 0.680000 1.0

pair_coeff 3 6 1.000000 1.0

pair_coeff 3 7 1.000000 1.0

pair_coeff 4 4 0.680000 1.0

pair_coeff 4 5 0.840000 1.0

pair_coeff 4 6 1.000000 1.0

pair_coeff 4 7 1.000000 1.0

pair_coeff 5 5 0.840000 1.0

pair_coeff 5 6 0.680000 1.0

pair_coeff 5 7 0.680000 1.0

pair_coeff 6 6 1.000000 1.0

pair_coeff 6 7 1.000000 1.0

pair_coeff 7 7 1.000000 1.0

bond_style harmonic

bond_coeff 1 27.6125 1.0

bond_coeff 2 27.6125 0.787234

angle_style cosine/squared

angle_coeff 1 2.5 180

angle_coeff 2 2.5 120

dielectric 20.0

special_bonds fene

neighbor 0.5 bin

neigh_modify delay 0 every 1 check yes

fix 1 all nve

fix 2 all momentum 1 linear 1 1 1 angular

fix 3 all press/berendsen iso 0.001250 0.001250 0.56045 modulus 4.465977

fix 5 oil temp/berendsen 0.498868 0.498868 0.56045

compute oil_temp oil temp

compute vpress all pressure thermo_temp virial

compute kepress all pressure thermo_temp ke

dump 1 all custom 2500 simulation1.movie id mol x y z

thermo 500

thermo_style custom step temp c_oil_temp pe press c_vpress vol

min_style cg

minimize 0.0 1.0e-8 1000 1000

velocity oil create 0.01 5334535 rot yes mom yes

timestep 0.011214

run 25000

##################Hold TEMP/PRESS (Stablize)##########################################

unfix 2

fix 2 all momentum 10 linear 1 1 1 angular

unfix 3

fix 3 all press/berendsen iso 0.001250 0.001250 0.56045 modulus 4.465977

unfix 5

fix 5 oil temp/berendsen 0.498868 0.498868 0.56045

dump 2 all custom 2500 simulation2.movie id mol x y z

thermo 500

thermo_style custom step temp c_oil_temp pe press c_vpress vol

timestep 0.011214

run 1000000

##################Hold TEMP/PRESS 160nsec (No Momentum Control)####################################

unfix 2

unfix 3

fix 3 all press/berendsen iso 0.001250 0.001250 0.56045 modulus 4.465977

dump 3 all custom 2500 simulation3.movie id mol x y z

thermo 500

thermo_style custom step temp c_oil_temp pe press c_vpress vol

timestep 0.011214

run 1000000

##################Analyze 80nsec #########################################################

dump 4 all custom 12500 run.txt id type x y z vx vy vz

dump 5 all custom 2500 simulation4.movie id mol x y z

compute msd_total all msd com yes

compute msd_oil oil msd com yes

thermo 1000

thermo_style custom step temp c_oil_temp c_msd_total[1] c_msd_total[2] c_msd_total[3] c_msd_total[4] c_msd_oil[1] c_msd_oil[2] c_msd_oil[3] c_msd_oil[4] pe press c_vpress vol

restart 500000000 simulation*.restart

timestep 0.011214

#ANALYZEDATA

run 500000000