read_dump with triclinic boundary conditions

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

Something seems to go wrong when using ‘read_dump’ to continue a simulation with triclinic boundary conditions in Lammps.

If the original dump file has a bounding box of the form

ITEM: BOX BOUNDS xy xz yz pp pp pp
-4.34705 10.9308 -4.14048
-0.15207 8.81237 1.72421
0.0669764 29.3269 -3.53651e-15

then if I dump another trajectory immediately after the ‘read_dump’ command with the ‘dump’ and ‘run 0’ commands I get

ITEM: BOX BOUNDS xy xz yz pp pp pp
-8.48753 12.655 -4.14048
-0.15207 8.81237 1.72421
0.0669764 29.3269 -3.53651e-15

Notice the change in the first two numbers. I would expect they should be left invariant under this operation.

According to the manual (http://lammps.sandia.gov/doc/dump.html) the format of the bounding box for triclinic boxes is

ITEM: BOX BOUNDS xy xz yz xx yy zz
xlo_bound xhi_bound xy
ylo_bound yhi_bound xz
zlo_bound zhi_bound yz

and the first two numbers are xlo_bound and xhi_bound. It seems that upon reading in the dump file in read_dump.cpp the numbers are being interpreted as xlo, xhi, etc, while I would think that they should be interpreted as xlo_bound, xhi_bound, etc. The relation between the two is given in http://lammps.sandia.gov/doc/Section_howto.html or domain.cpp.

In the VMD routine that processes dump files (lammpsplugin.c) a problem exists as well (already communicated to Axel), as the conversion between xlo_bound and xlo is slightly different than the one given in the manual.

Thank you,
Bart

2

the recommended way to continue a simulation is to use
write_restart and read_restart.

axel.

3

Fair point, although a warning and/or error message might then be
appropriate.

The reason why I was actually using read_dump was that restart2data was
giving a segmentation fault.

A closer look shows that the fault can be circumvented by simply
commenting out line 447 that calls autodetect in restart2data.cpp. The
autodetect function opens the 'restartfile' instead of the 'basefile'
variable in restart2data.cpp when using a file with a % sign in it.

Thanks,
Bart

4

Well,
It should still work correctly.
if you want to have some fun you could use the user-molfile package and the VMD plugin to read the coordinates.

Axel

5

Think I fixed this issue in restart2data - will be in next patch.

I'll look at read_dump for triclinic
boxes - we may be doing something wrong there.

Thanks,
Steve

6

Thank you both for the replies, and for the fix. That package Axel
mentioned might be useful for me, I'll have a closer look at it.

Sorry to bug you further, but it seems that there are some more problems
with triclinic boundary conditions. I guess it is a feature not used very
often.

I encountered it by trying to set the center of mass (COM) of the system
to zero by:
"
variable ax equal -xcm(all,x)
variable ay equal -xcm(all,y)
variable az equal -xcm(all,z)
print "-com= \{ax\} {ay} \{az\}" displace\_atoms all move {ax} \{ay\} {az} units box
variable ax equal xcm(all,x)
variable ay equal xcm(all,y)
variable az equal xcm(all,z)
print "com \{ax\} {ay} ${az}"
"

One would expect that the last print gives three numbers very close to
zero. However, a test run showed values substantial deviating from zero
for the x and y component.

The problem seems that the unboxing is not properly accounted for. For
example, in group.cpp function xcm has the code

cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
        cmone[2] += (x[i][2] + zbox*zprd) * massone;

which gives the non-zero COM values for triclinic boundary conditions with
non-zero tilt. Changing it to

cmone[0] += (x[i][0] + xbox*xprd+ybox*xy+zbox*xz) * massone;
        cmone[1] += (x[i][1] + ybox*yprd+zbox*yz) * massone;
        cmone[2] += (x[i][2] + zbox*zprd) * massone;

resolves it. Note that a similar correct implementation has been applied
in, e.g., void Domain::unmap(double *x, tagint image) of domain.cpp, void
ComputeDisplaceAtom::compute_peratom(), compute_displace_atom.cpp, and
void FixRigid::setup(int vflag) in RIGID/fix_rigid.cpp

group.cpp is not the only file that contains code that appears to be
affected. Also the following files might need further inspection, as they
use a similar construct:
- compute_com_molecule.cpp
- compute_gyration.cpp
- compute_gyration_molecule.cpp
- compute_msd_molecule.cpp
- fix_addforce.cpp
- fix_momentum.cpp
- fix_spring_rg.cpp
- fix_spring_self.cpp
- fix_tmd.cpp
- velocity.cpp
- POEMS/fix_poems.cpp
- USER-MISC/fix_addtorque.cpp
- USER-MISC/compute_temp_rotate.cpp

Thank you,
Bart

7

Thank you both for the replies, and for the fix. That package Axel
mentioned might be useful for me, I'll have a closer look at it.

Sorry to bug you further, but it seems that there are some more problems
with triclinic boundary conditions. I guess it is a feature not used very
often.

i don't think it is so much the frequency of use, but rather
the fact that non-orthogonal cells is a feature added at a
later point and that new features are usually only checked
against old features that the person implementing it needs.

we don't have a culture of systematic unit and regression
testing and thus a lot of changes can lead to problems
that can go unnoticed for extended periods of time, if a
certain combination of features is not used by somebody
that is as careful in checking results as you are (most
people are not). as a person maintaining a few 'user'
packages, i was 'bitten' by this kind of attitude several
times already and whenever i have the time, i am carefully
reading through the commit diffs to lammps repository in
order to make sure, i won't miss anything.

that being said, testing for all combinations of features
in a software as flexible and capable as lammps is near
impossible. quite a few issues, for example, only arise
when you run on parallel, or have an unusual distribution
of particles and so on.

having somebody going over things in such a systematic
way, is a very helpful and commendable thing.

cheers,
     axel.

8

Fixed this with the 23 Oct patch today - it was indeed a bug.

Thanks,
Steve

9

Think I just changed all of these to work
with triclinic geometries. Will post a patch,
thanks.

Steve