[lammps-users] Velocity Verlet algorithm wrong? // change of angle parameters at runtime

Dear all

I have started to work with LAMMPS 3 months ago and have studied the source
code since then. Right now i got two questions:

1. The first one is probably quiet easy to answer. I guess i didnt
understand some part of the code or overlooked something critical. My
understanding is that LAMMPS uses a velocity verlet algorithm if I do a run
using the fix_NVE command. If i break up the algorithm as I learned it, it
should do this in one timestep:

  1) update v (velocity first half step)
  2) update x (locations)
  3) update v (velocity second half step)
  4) update a (forces or acceleration)

(see this page for example:

Studying the Lammps source code I found that it did something else,
something I dont really understand.

The main loop is the verlet::iterate routine. The first point something is
updated is the modify->initial_integrate() command.
Following this routine I find it just calls the initial_integrate of the fix
I use. In the fixnve::inital_integrate I find it updates first the
velocitys, and then the positions.

Back to the verlet::iterate routine one sees that the next thing that is
updated are the forces (force->pair->compute or whatever forces one uses).
After these are updated, the modify::final_integrate routine is called,
which in turn calls the final_integrate of the fix one uses. Here the
velocities are updated again.

To summarise it it looks for me as if LAMMPS does the following order:
  1) update v
  2) update x
  3) update a
  4) update v

This is not the same order as I learned it for the velocity verlet
algorithm. And I dont see that step 4 would be necessary as the code is now.
If I just put a factor 2 in the update formula in the initial_integrate
routine and do not use the final integrate routine, Lammps should do exactly
the same as it does now. Any clarification would be appreciated much.

2. For a Li-B-O system I need to use angle potentials. Unfortunately the
coordination of the B atoms may change during the simulation. Therefor I
need to change how many angles each B atom has, which atoms are associated
to which angle and what parameters the angle potentials of a certain B atom
have (the parameters are different for a triple or quadruple coordination).
Did anyone do something like that allready?
I started to work on that problem this week. My idea would be to give every
B atom 6 angles right from the start and to use 3 different angle types. One
for the triple coordination one for the quadruple coordination and one as a
place holder with parameters equal zero. Now i would start messing around
with the neighbor->anglelist, atom->angle_type and atom->angle_atom1/2/3
lists at runtime. Updating these lists with the situation I find every n
Any comments on that idea? Have i forgotten any lists I would need to

Thank you for your help, and please forgive this long text right in my first

Christian M�ller

P.S. I implemented allready a routine that calculates density distributions
for the different atom types over the simulation time. If someone else needs
that feel free to mail me.

  1) update v (velocity first half step)
  2) update x (locations)
  3) update v (velocity second half step)
  4) update a (forces or acceleration)

This velocity-Verlet is not correct, b/c you can't do the 2nd v update
until you have new forces. So you need to flip steps 3/4 which is
what LAMMPS does and what the algorithm on the WWW page
you mention does. I.e the WWW page matches LAMMPS, so I don't
follow why you think your version matches the WWW page.

Also, in practical terms you can't delay doing step 4 until
the beginning of the next timestep, else you will not
have a consistent set of x,v at the end of the timestep to
do analysis on or dump out.

Re: your Li-B-O system with changing angles. You have a couple
choices. (1) You could implement it as a pair potential (no angles)
that looks for 3-body triplets each timestep and computes the
angle terms. This is what bond-order potentials like Tersoff (for Si)
do. (2) You could break and form angles in the angle potential.
But this is quite tricky to do since you must update the LAMMPS data
structures you mention, and you must do it correctly in parallel (for
triplets that straddle processor boundaries), and you must do it
on the correct timesteps when data is about to be passed between
procs due to re-neighboring. So I advise against it. (3) Define
all possible angles and have the angle potential only compute the
currently valid ones. This would only work if your system is
static (a solid?) so that no new angles possibly form over time.


lammps users,

I have been writing a fix for lammps-30Nov05 and found an undocumented
change (I didn't see it mentioned on the web page) in lammps-12Apr06
for a method in the Fix class virtual method memory_usage():

30Nov05 - Fix::memory_usage(int)

12Apr06 - Fix::memory_usage(void)

Without making the change, lammps will compile but hang during
execution because a memory_usage(int) method is syntactically acceptable,
but lammps internals are expecting a memory_usage(void) method, and
the default is to return 0;

I hope this brief note is helpful.


Yes - the call to memory_usage() in fixes did undergo a syntax
change. This was part of the 15 Mar 06 feature note, though
the notice on the bug page didn't talk about this level of detail.

When I compile a fix with the old memory_usage(void),
I get a compiler warning that I am overridiing the new memory_usage(int)
base class version, so that should indicate the problem ...