Problem with lipids all moving in same direction during NPT

LAMMPS LAMMPS Mailing List Mirror
1

Hello all,

I am attempting to equilibrate, and eventually deform, a POPC lipid bilayer. However, I am running into a problem where all of my lipid molecules begin moving in the same direction during NPT. This movement is amplified when I attempt to deform the bilayer. This appears to be independent of lipid bilayer size and water amount.

I am using the 9Aug13 version of LAMMPS. I am currently running on 96 processors, but previously used 192 processors.

I do not understand where the motion could be coming from. Am I missing something obvious that would cause the molecules to begin moving?

I have attached an animated gif displaying the top of the lipid bilayer. The water molecules are hidden. The lipid head group of one phospholipid has been changed to a different color to aid in visualization of movement.

Additional Information that may be useful:

Initial Structure:
I retrieved a 72 lipid (36/leaflet) POPC bilayer structure in pdb format from the site http://www.lobos.nih.gov/mbs/index.shtml which I converted to a LAMMPS data file using the charmm2lammps tool, a psf file produced in VMD, and the forcefield parameters for CHARMM36. I then changed the parameter values as directed in the LAMMPS documentation to use TIP3P water with long-range Coulombic solver and added more water. The structure has 36,858 atoms with the extra water.

Some Simulation Info:
I begin by initializing the velocities, moving the x and y dimensions of the box in to more closely fit the sides of the lipid bilayer, and ramping the temperature to 310 K while running NVT. I then change to NPT and run the simulation for the remainder of 3 ns (simulation is then restarted to run another 3 ns). I am using a .5 fs timestep. I have included portions of my input file:

---------- Define Simulation Parameters ------------------

units real
neigh_modify delay 0 every 1 check yes

atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic

pair_style lj/charmm/coul/long 8 10
pair_modify mix arithmetic
kspace_style pppm 1e-4
special_bonds charmm

---------- Define Thermo Style ------------------

thermo 10
thermo_style custom step elapsed cpu temp etotal pe ke press evdwl ecoul ebond eangle edihed eimp elong vol fmax fnorm

---------- NPT Portions of Simulation ------------------

#Uncoupled box with pressure = 1 atm and temperature = 310 K
fix 5 all npt temp 310 310 100 aniso 1 1 500

#Use “every” command without options during NPT to force kspace grid to reset
run 5700000 every 1000 “”

Results:
The Volume, Pressures, Energies all appear to become equilibrated (no net increase/decrease) during NPT. Additionally, the area per lipid appears to be acceptable. However, all of the lipids develop a motion in the same direction.

Other Comments:
Under advice of one of my committee members, I also ran a short NVE simulation starting after one of my equilibration runs (after 6 ns) which resulted in the potential, kinetic, and total energies having a slow increase.

I can include the complete input file or a data file (~8MB) if needed.

Thanks,
Michael

Top_no-water.gif

2

Hello all,

I am attempting to equilibrate, and eventually deform, a POPC lipid bilayer.
However, I am running into a problem where all of my lipid molecules begin
moving in the same direction during NPT. This movement is amplified when I
attempt to deform the bilayer. This appears to be independent of lipid
bilayer size and water amount.

looks a lot like you run into the early stages of what is called the
"flying icecube syndrome".

during equilbration your system could be picking up some center of
mass motion and that does not easily dissipate without a dissipative
thermostat (which npt is not).

there are a number of contributing factors:
- small cutoffs
- not tight enough convergence on either pppm or shake or both
- too large time step
- modifications to the system that result in inducing a center of mass momentum.

since you already report some lack of energy conservation, you should
work on improving that, but also you can use the velocity command
after the nvt run to remove any center of mass momentum present at
that time. it may also help to translate the system, so that the
origin becomes the center of the box, if that is not already the case.
the sum of these steps should improve the situation and help you to
identify, how deformation will impact the center of mass motion.

that, if it is still present, may then be quelled through means like
using fix momentum. be careful to only use it after you took care of
all other options, since it may otherwise mask residual problems.

axel.

3

Hey Axel,

Thanks for the comments.

I will try some simulations based on your suggestions and try to respond again in a few days.

I apologize if this is an obvious question, but how would I determine if the pppm convergence is tight enough?

Michael

4

Hey Axel,

Thanks for the comments.

I will try some simulations based on your suggestions and try to respond
again in a few days.

I apologize if this is an obvious question, but how would I determine if the
pppm convergence is tight enough?

well, what is "good enough" in general?

a convergence of 1.e-4 is rather loose and sensitivity to errors is
larger for variable cell calculations.

for testing, you can compare forces to a very tightly converged value
and also compare total energies.

axel.

5

Hey Axel,

Resetting the linear momentum with the velocity command after the NVT portion of the simulation seems to have fixed the problem.

I am still planning to run an NVE simulation to be sure of energy conservation.

Thanks for the suggestions,

Michael Murphy

6

Hey Axel,

I managed to maintain energy conservation after some modifications to the kspace settings (kspace_style pppm 1e-4 with kspace_modify ad).

However, after I began performing equilibrations again, I have noticed that the motion begins to develop again ~1 ns after setting the linear momentum to zero while running NPT. The new motion development only occurs in some of the simulations. Is it likely these motions are a result of numerical errors accumulating? Or is it more likely some other aspect in my simulation is initiating the motion and needs to be corrected?

I am still using the 9Aug13 version of LAMMPS on 96 processors.

Thanks,

Michael Murphy

7

Hey Axel,

I managed to maintain energy conservation after some modifications to the kspace settings (kspace_style pppm 1e-5 with kspace_modify ad).

However, after I began performing equilibrations again, I have noticed that the motion begins to develop again ~1 ns after setting the linear momentum to zero while running NPT. The new motion development only occurs in some of the simulations. Is it likely these motions are a result of numerical errors accumulating? Or is it more likely some other aspect in my simulation is initiating the motion and needs to be corrected?

I am still using the 9Aug13 version of LAMMPS on 96 processors.

Sorry for the second email. Noticed I didn’t update the kspace value in the first.

Thanks,

Michael Murphy

8

Hey Axel,

I managed to maintain energy conservation after some modifications to the
kspace settings (kspace_style pppm 1e-5 with kspace_modify ad).

However, after I began performing equilibrations again, I have noticed
that the motion begins to develop again ~1 ns after setting the linear
momentum to zero while running NPT. The new motion development only occurs
in some of the simulations. Is it likely these motions are a result of
numerical errors accumulating? Or is it more likely some other aspect in
my simulation is initiating the motion and needs to be corrected?

if your simulation keeps a zero total momentum for about 1ns that is
already a success. at this level some error accumulation is to be expected.
you may reduce it by choosing the origin and by perhaps using an even
tighter pppm convergence. the ad method seems to be a bit less accurate at
the same convergence parameter.

you should be able to quell the remaining COM drift using fix momentum..

axel.