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:
http://www.ph.ed.ac.uk/~graeme/compmeth/verlet.html)
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
timesteps.
Any comments on that idea? Have i forgotten any lists I would need to
update?
Thank you for your help, and please forgive this long text right in my first
mail.
Yours
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.