[lammps-users] Pressure difference in subsequent runs

Hello,

I noticed that when using two "run" commands right after one another, the
pressure reported at the end of the first run is different from the
pressure at the beginning of the second run. For example, when I run the
attached input file, I get the following output (suitably cropped):

Step Temp E_pair E_mol TotEng Press
     500 1.0142374 70.327094 0 71.847707 23.948199

Step Temp E_pair E_mol TotEng Press
     500 1.0142374 70.327094 0 71.847707 23.371509

I usually split up long simulation runs into several shorter ones to obtain
the simulation statistics that is written after each run command, so I can
check for dangerous neighbor list builds before the long simulation is
completed. This seemingly small difference in pressure is very noticable
in large systems.

I also have another (and possibly related) question: in DPD simulations, is
the velocity-dependend drag force included in the pressure calculation?

Thank you very much,

  Lutz

dpd.in (521 Bytes)

We have tested the following input scripts for LAMMPS but we couldn't
perform an NPT simulation for a crystal with a surface. In fact we are
trying to perform a thermally expansion for the crystal with surface.
How can we simulate the expansion for a vacuum with "kspace_style" and
"kspace_modify" commands? Since NPT needs periodic boundaries to control
the pressure, the ewald with slab option couldn't initialized. Thanks
for any help.

Best regards.

  Berk ONAT
  Informatics Institute
  Istanbul Technical University

input script

Berk,

Try turning off the periodicity and pressure control in the z-dimension and it should work.

See:
http://lammps.sandia.gov/doc/kspace_modify.html
http://lammps.sandia.gov/doc/boundary.html
http://lammps.sandia.gov/doc/fix_npt.html

Something like:

boundary p p f

and:

fix 1 all npt 50.0 50.0 0.04 aniso 0.0 0.0 0.0 0.0 NULL NULL 0.1

Also, you may want to use PPPM instead of Ewald for better speed.

Paul

Lutz,

To answer your last question first: yes, the velocity-dependend drag force is indeed included in the pressure calculation. (See lines 125 - 163 of pair_dpd.cpp, especially lines 126 and 130.)

So your questions are related. Since dpd includes a random force (see lines 127 and 131), and the virial includes all of the forces, it is not surprising that the pressure comes out different each time it is calculated, even when the velocities and positions of the atoms don’t change.

Paul

Paul,

thank you very much for your reply. I think that the difference in the
reported pressure between subsequent runs is not only due to different
realizations of the random forces. I attach the thermo output of a series of
40 runs, each consisting of 2000 timesteps.

If you plot the pressure as a function of time (for example, 'plot
"log.lammps" u 1:6 w lp' in gnuplot), you see spikes where the pressure dips
below its average every 2000 timesteps. These spikes correspond to the first
thermo output of each run. Maybe the random forces have not been calculated
yet at this stage of a run?

  Lutz

log.lammps (70.8 KB)

So many thanks Paul. We have changed the script and it works for
"boundary p p f" now. But we still have a problem. We can not use EAM
potential with kspace_style Ewald or PPPM. I have double checked LAMMPS's
documents and sourceforge site, but couldn't find any information about
this problem. Is there any way to use EAM potentials with PPPM or Ewald?

  Best regards.

  Berk

EAM doesn't store or use charge on atoms, so what does
it mean to use it with PPPM ?? In other words,
EAM is a complete potential, with parameters fit to a cutoff;
there is no long-range Coulombic component.

Steve