Freezing and Drifting Problem

LAMMPS LAMMPS Mailing List Mirror
1

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

2

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.

as far as i recall, this is a known issue of the martini model.
for water, it requires special "anti-freeze" particles and also
requires the use of a fairly large time step. the shifted freezing
point of water is common to all similar coarse grain models.
with 4 waters per particle it more severe than with 3 waters
per particles and so on. you have to consider your system
as "partially pre-frozen".

i have no other recommendation than switching to a different
coarse grain model that suffers less from this problem, but
then again, i am a bit biased. :wink:

cheers,
     axel.

(stupid question: why on earth did you switch from the much
more convenient "real" units to reduced units??)

3

I have heard that there is a freezing problem, but I had not heard of the drifting problem.
I chose the martini model haphazardly. I really just want to complete this problem and move on.
I went non-dimensional a while ago in order to better understand the system and compare it to other models.
The time step you see is equivalent to 20 fs.

4

I have heard that there is a freezing problem, but I had not heard of the drifting problem.

those are connected. the kinetic energy has to go somewhere.
you may have heard the term "flying icecube syndrome".

using a berendsen thermostat makes it worse, since its
adjustment of the kinetic energy through global velocity scaling
does effectively enhance the drift. using fix langevin with
the zero yes flag set may be a better choice to counter
any potential drift.

axel.