Kunwoo,
(1)Velocity setting at initial stage
I ran a simulation without setting "velocity_create"
for npt simulation and the average temperature was
323K as I set using "fix_npt"(in.HD-n). My first
question is if the "velocity_create" command affects
the final result of simulation. I think as long as
the temperature is maintained constant, it doesn't
affect the average thermodynamic properties. But I
just want to make sure if my thought is right.
See the LAMMPS documentation here:
http://www.cs.sandia.gov/~sjplimp/lammps/doc/Section_errors.html#9_1
Pasted for your convenience here:
"If two LAMMPS runs do not produce the same answer on
different machines or different numbers of processors,
this is typically not a bug. In theory you should get
identical answers on any number of processors and on
any machine. In practice, numerical round-off can
cause slight differences and eventual divergence of
molecular dynamics phase space trajectories within a
few 100s or few 1000s of timesteps. However, the
statistical properties of the two runs (e.g. average
energy or temperature) should still be the same.
If the velocity command is used to set initial atom
velocities, a particular atom can be assigned a
different velocity when the problem on different
machines. Obviously, this means the phase space
trajectories of the two simulations will rapidly
diverge. See the discussion of the loop option in the
velocity command for details."
So, yes, in theory it shouldn't matter what your
initial conditions are as long as you allow the system
to equilibrate, are sampling at steady-state in the
correct ensemble, and perform adequate sampling. In
this sense, the statistical (thermodynamic) properties
should be independent of the initial conditions.
(2)Initial temperature
My another question is related to the initial
temperature. I found that the two runs with or
without "velocity_create" resulted in the same
initial temperature(step 0) of 416K which is not the
one I set in my input file(in.HD-d and in.HD-n). I
don't know why did that happen. Does anyone know
what is wrong with my input file?
If you're using shake, degrees of freedom are removed
from the system due to the shake constraints. This in
turn affects the temperature calculation. The
velocity_create command assumes there are no
constraints in the system when it assigns velocities
to get the user-specified temperature. Again, this
shouldn't matter much as long as you're using a
thermostat and give the system a bit of time to
equilibrate. If you assign velocities in your data
file that happen to yield the same initial T as you
get with the velocity_create command, that would
explain why you see the same T irregardless of whether
or not you use the velocity_create command.
(3)fix_npt
I ran two different simulations with the same input
file except the fix_npt command(and the
velocity_create command too). In one run, I set the
lateral pressure(Pxx=Pyy) 1.0 and the normal
pressure(Pzz) 1.0 separately using "fix_npt xy" as
shown in in.HD-d file. On the other run, I just set
the external pressure 1.0 in all three directions
using "fix_npt xyz" command as shown in in.HD-n
file. My results show that the thermodynamic
properties are not quite different, but I found
noticeable difference in some structural difference,
for example the run using "fix_npt xy" resulted in
~30% increase in the lipid head group area. So, is
it quite common or something wrong with my results?
There is a large and complex literature about whether
to run lipid bilayer simulations with constant surface
area or constant surface tension. In our work (see
reference below), we've opted for the simpler and more
straight-forward constant membrane surface area
approach, which I would recommend to you. Using the
number of lipids you have in each leaflet and the
experimental surface area per lipid, you compute the
desired membrane surface area and keep it constant
throughout the simulation. I'd recommend setting the
lateral dimension pressure to 1 atm and allowing your
simulation cell to fluctuate in that dimension. Here's
a reference to a paper we've written that refers to
some of the literature on this topic:
Crozier, P.S., Stevens, M.J., Forrest, L.R., Woolf,
T.B., �Molecular dynamics simulation of dark-adapted
rhodopsin in an explicit membrane bilayer: Coupling
between local retinal and larger scale conformational
change,� J. Mol. Biol., 333, 493-514, (2003).
See below for some settings that I typically use
(especially note the fix npt settings).
Paul
units real
neigh_modify delay 4 every 1
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
pair_style lj/charmm/coul/long 8 10
pair_modify mix arithmetic
kspace_style pppm 1e-4
read_data m51.data
special_bonds charmm
fix 1 all shake 1e-6 500 0 m 1.0 a 36
fix 2 all npt 298. 298. 100.0 aniso NULL
NULL NULL NULL 1.0 1.0 1000.0 drag 1.0
thermo 500
timestep 2.0
restart 50000 51.res
dump 1 all atom 500 m51.dump0
run 500000