[lammps-users] Problems running basic minimize-heat-equilibriate simulation

LAMMPS LAMMPS Mailing List Mirror
1

Hi,
I'm a new LAMMPS user. My group at Wayne State University in Detroit, Mich. researches phospholipid bilayers, among other things and has largely used NAMD for past modeling.

I've been trying to get LAMMPS running with some basic simulations to prove out that its working as intended, before moving on to more complex simulations (phospholipids).

I have a working simulation in NAMD that simulates 300 ethane molecules at the saturation pressure for 220 K (appr. 4.9 atm). The goal of the simulation is to get an equilibriated density, and make sure it matches then NIST experimental data. This simulation involves three configuration files -- a minimization one, a heating one, and finally an equilibriation one. Having achieved good results in NAMD -- a equilibrium density of 16.29 mol/l (experimentally predicted is 16.5), I'm now trying to translate this into LAMMPS.

I ran the charmm2lammps perl script to extract my .data file from my PDB/PSF model, and then modified the input file (see attached for these files). I had to add a minimization step (expected) or otherwise the simulation would fail as atoms would go out of range, likely due to explosion. That brings me to my first question -- what kind of minimization values are generally a good idea to use (my system has 300 molecules, 2400 atoms total)??

Using a minimization, I next tried to take the simulation through a heating process at constant pressure, using the nph and temp/rescale fixes. This did not work, the simulation exploded.

I then tried directly jumping into my desired simulation temperature 220 K via the velocity command (previously I had been calling this command with a temp of 0 K, I now am using 230 K). I then run a million step simulation with a npt fix with the saturation pressure and temp of 220 K. My dampening constants for temp and pressure, respectively, are 100 and 1000 fs.

This simulation runs, but it does not go to the liquid phase, even after a long minimization. I packed the system to approximately the appropriate volume as well -- 27,000 Angstrom^3. I've tried turning up the pressure and running at above the critical temperature, but no matter what the system seems to explode as the Potential energy and volume soar up.

Looking at my files, does anyone have an idea how I can achieve my object -- to simulate this box of ethane in the liquid phase @ 220 K? What am I doing wrong?

Also, if anyone knows of good LAMMPS tutorials for beginners that they can point me to, that would be great. The limited examples that came with the distribution of LAMMPS weren't overly helpful. I believe I understand most of the commands, but given my previous difficulties, putting them together into a working simulation still seems tedious.

I appreciate any help, thanks!!

Cheers,
Jason Mick
Graduate/Ph.D Student
Wayne State University
E: [email protected]...
C: 248.978.5941

ethane_30A3_3.in (757 Bytes)

ethane_30A3_3.pdb (185 KB)

2

I don't see anything obviously wrong with what you've done. It would
be good to check that the CHARMM translation you've done is
giving an equivalent force field. You should be able to
compare energies for one or more snapshots in both LAMMPS and NAMD and
get essentially identical numbers. If things are blowing up, then
I'd look at the log file and viz (doing thermo output frequently) to see
what is going bad - possibly a big volume fluctuation or the like.

Paul - do you have any other ideas on this?

Steve

3

Few questions.

  1. Did you cross-check that the potential parameters that lammps is using are same as you want them to be?
  2. How do you know that it is not liquid at 220 K?

It will be good see the temperature, pressure and volume evaluation of your system to see whether all of them have got equilibrated after one million steps. Few graphs may be appropriate for the case of npt simulations at 220 K and assigned 4.9 pressure. Also, do you really want to save restart file every 10 steps?

Vikas

4

Jason,

I agree with what Steve and Axel have said on this. Specifically, I'd second Axel's suggestion to check the box size. I think that the charmm2lammps script may not yield exactly the right box size that you had in NAMD, so you might have to change that by hand to get agreement. Once you've checked the box size and force field parameters, I agree with Steve that checking that both NAMD and LAMMPS give the same timestep 0 energies is a really good sanity check. Given the same timestep size and simulation parameters, subsequent timesteps should also produce the same energies, but gradually diverge due to small numerical differences between the codes.

Paul

5

parameters, I agree with Steve that checking that both NAMD and LAMMPS
give the same timestep 0 energies is a really good sanity check. Given
the same timestep size and simulation parameters, subsequent timesteps
should also produce the same energies, but gradually diverge due to
small numerical differences between the codes.

the differences in trajectories between NAMD and LAMMPS are more
likely to be larger than smaller for two reasons: the random number
generation in the langevin thermostat and - when running in parallel -
the load balancing and different division into subdomains (patches
in "NAMD-speak") between the two codes due to the non-associativity
of regular floating point math.

cheers,
   axel.