Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix for issue #3383 triclinic lammpsdump #3403

Merged
merged 9 commits into from
Mar 23, 2022
Merged

Fix for issue #3383 triclinic lammpsdump #3403

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
5c10286
changes DumpReader to expect orthogonoal bounding boxes of triclinic …
alexgorfer Aug 24, 2021
e7c5af4
test between data, dump and expectation
alexgorfer Dec 19, 2021
ad456c7
pep8 fixes
alexgorfer Dec 19, 2021
ef2ef22
comments and assert_allcose
alexgorfer Mar 8, 2022
c2c331c
Merge branch 'MDAnalysis:develop' into origin/issue-3383
alexgorfer Mar 9, 2022
3dd0cdd
doc and changelog
alexgorfer Mar 10, 2022
5b14120
Merge branch 'develop' into origin/issue-3383
alexgorfer Mar 20, 2022
deb9e64
Changelog, IssueLink, pep8
alexgorfer Mar 23, 2022
5056fd5
pep8
alexgorfer Mar 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions package/MDAnalysis/coordinates/LAMMPS.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -560,10 +560,16 @@ def _read_next_timestep(self):

triclinic = len(f.readline().split()) == 9 # ITEM BOX BOUNDS
if triclinic:
xlo, xhi, xy = map(float, f.readline().split())
ylo, yhi, xz = map(float, f.readline().split())
xlo_bound, xhi_bound, xy = map(float, f.readline().split())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be possible to add a relevant test (maybe based on the issue this aims to fix) that covers these lines?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah of course, I can write a test that reads in a triclinic box. This would test that it doesnt fail. Testing that the actual issue would be resolved is difficult/impossible?. I think one would have to create a triclinic box in lammps, then actually run lammps and read the output to see if the dimensions are the same.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So here that's what we would require to do - create a triclinic box in lammps and then check that we reproduce expected outcomes (calculated by hand and/or comparing to a lammps output).

If we can't do such a correctness check then we might have to discuss further whether or not we can offer this option.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have a first correctness check here:
triclinic_lammps_2021_09_02.zip
I use the values from the input data file for lammps that require xlo, xhi etc. directly and compare those dimensions to the ones imported from the lammpsdump file of the same system:

import MDAnalysis as mda
import numpy as np
from MDAnalysis.lib import mdamath

u = mda.Universe("dump.atom", topology_format='LAMMPSDUMP')

## Direct values from lammps data file
xlo = -0.32115478301032807
xhi = 16.831069399898624
ylo = -0.12372358703610897 
yhi = 25.95896427399614
zlo = -0.045447071698045266 
zhi = 12.993982724334792
xy = 1.506743915478767
xz = - 6.266414551929444 
yz = - 0.42179319547892025

box = np.zeros((3, 3), dtype=np.float64)
box[0] = xhi - xlo, 0.0, 0.0
box[1] = xy, yhi - ylo, 0.0
box[2] = xz, yz, zhi - zlo

direct_dimensions = mdamath.triclinic_box(*box)

print(direct_dimensions == u.dimensions)

resulting in the same dimensions

[ True  True  True  True  True  True]

Would this be sufficient? Should I then write this procedure into a test file or what would be the best practice for testing specifications set in other software?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can write something like this into the testsuite! the relevant file is here.

You will need to learn a bit of pytest, and use the text fixtures but hopefully the file is commented enough to be comprehensible. The numpy testing functions will be very helpful. See here for how to run the testsuite.

If you need to add an additional file to the test data, please do add it to the mdanalysis/testsuite/MDAnalysisTests/data/lammps folder and import it in
mdanalysis/testsuite/MDAnalysisTests/datafiles.py

I am going to tag @richardjgowers in case I am missing something about ortho vs triclinic in MDA/LAMMPS here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah awesome thank you! I will get the keywords you recommended and then write it up in there!

hmacdope marked this conversation as resolved.
Show resolved Hide resolved
ylo_bound, yhi_bound, xz = map(float, f.readline().split())
zlo, zhi, yz = map(float, f.readline().split())

# converts orthogonal bounding box to the conventional format, see https://docs.lammps.org/Howto_triclinic.html
hmacdope marked this conversation as resolved.
Show resolved Hide resolved
xlo = xlo_bound - min(0.0, xy, xz, xy + xz)
xhi = xhi_bound - max(0.0, xy, xz, xy + xz)
ylo = ylo_bound - min(0.0, yz)
yhi = yhi_bound - max(0.0, yz)
hmacdope marked this conversation as resolved.
Show resolved Hide resolved

box = np.zeros((3, 3), dtype=np.float64)
box[0] = xhi - xlo, 0.0, 0.0
box[1] = xy, yhi - ylo, 0.0
Expand Down
48 changes: 45 additions & 3 deletions testsuite/MDAnalysisTests/coordinates/test_lammps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@
import MDAnalysis as mda
from MDAnalysis import NoDataError

from numpy.testing import (assert_equal, assert_almost_equal)
from numpy.testing import (assert_equal, assert_almost_equal, assert_allclose)

from MDAnalysisTests import make_Universe
from MDAnalysisTests.coordinates.reference import (
RefLAMMPSData, RefLAMMPSDataMini, RefLAMMPSDataDCD,
)
from MDAnalysisTests.datafiles import (
LAMMPScnt, LAMMPShyd, LAMMPSdata, LAMMPSdata_mini, LAMMPSDUMP,
LAMMPSDUMP_allcoords, LAMMPSDUMP_nocoords
LAMMPScnt, LAMMPShyd, LAMMPSdata, LAMMPSdata_mini, LAMMPSdata_triclinic,
LAMMPSDUMP, LAMMPSDUMP_allcoords, LAMMPSDUMP_nocoords, LAMMPSDUMP_triclinic
)


Expand Down Expand Up @@ -625,3 +625,45 @@ def test_auto_is_unscaled_match(self, universes):
for ts_a, ts_s in zip(universes["auto"].trajectory,
universes["unscaled"].trajectory):
assert_almost_equal(ts_a.positions, ts_s.positions, decimal=5)


class TestLammpsTriclinic(object):
@pytest.fixture()
def u_dump(self):
return mda.Universe(LAMMPSDUMP_triclinic, format='LAMMPSDUMP',
lammps_coordinate_convention="auto")

@pytest.fixture()
def u_data(self):
return mda.Universe(LAMMPSdata_triclinic, format='data',
atom_style='id type x y z')

@pytest.fixture()
def reference_box(self):
# manually copied from triclinic data file
xlo = -0.32115478301032807
xhi = 16.831069399898624
ylo = -0.12372358703610897
yhi = 25.95896427399614
zlo = -0.045447071698045266
zhi = 12.993982724334792
xy = 1.506743915478767
xz = -6.266414551929444
yz = -0.42179319547892025

box = np.zeros((3, 3), dtype=np.float64)
box[0] = xhi - xlo, 0.0, 0.0
box[1] = xy, yhi - ylo, 0.0
box[2] = xz, yz, zhi - zlo

return mda.lib.mdamath.triclinic_box(*box)

def test_box(self, u_dump, u_data, reference_box):
hmacdope marked this conversation as resolved.
Show resolved Hide resolved
# NOTE threefold validation testing both data and dump reader.
for ts in u_dump.trajectory:
assert_allclose(ts.dimensions, reference_box, rtol=1e-5, atol=0)

for ts in u_dump.trajectory:
assert_allclose(ts.dimensions, u_data.dimensions, rtol=1e-5, atol=0)

assert_allclose(u_data.dimensions, reference_box, rtol=1e-5, atol=0)
34 changes: 34 additions & 0 deletions testsuite/MDAnalysisTests/data/lammps/albite_triclinic.data
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
LAMMPS triclinic data file

17 atoms

1 atom types

-0.32115478301032807 16.831069399898624 xlo xhi
-0.12372358703610897 25.95896427399614 ylo yhi
-0.045447071698045266 12.993982724334792 zlo zhi
1.506743915478767 -6.266414551929444 -0.42179319547892025 xy xz yz

Masses

1 26.9815

Atoms # atomic

192 1 2.939929226745528 0.28126611328982504 0.509212291451447 0 0 0
85 1 0.2851050832641419 3.0807154102734917 4.6247763608193155 0 0 0
295 1 -0.8077822604729055 3.9571902907303733 3.9345748176468693 0 0 0
300 1 1.8093279024899684 3.3208268390257514 4.215038876613601 0 0 0
188 1 5.449513973050337 0.4263205981754459 1.587561248258666 0 0 0
191 1 3.999161625604918 5.667029705285251 1.8280363439039549 0 0 0
299 1 -0.12308953374879046 1.5325054443344792 4.1337384269477715 0 0 0
159 1 1.4500667066314719 1.1149430067523804 2.391995904640104 1 0 1
136 1 6.6253293198545045 2.5456109021590163 2.847256975857115 0 0 0
146 1 5.493605209601826 0.261491883863461 4.532083302884116 0 0 0
193 1 3.5376239452534324 3.671398031996834 -0.02020726596457809 0 0 0
81 1 -1.0940910841830676 5.6783700063066815 3.9223430559877266 0 0 0
189 1 5.995616400934193 4.020967239135711 0.7068582903868794 0 0 0
43 1 6.847965492945945 0.4349078018589978 0.7454921986075675 0 0 0
304 1 4.380319548944577 2.5314626579385897 5.491572411560825 0 0 0
86 1 -0.009824994256196633 0.1211322463682511 5.194921290704171 0 0 0
302 1 3.8457424411075216 5.270216285190885 4.93718329708941 0 0 0
26 changes: 26 additions & 0 deletions testsuite/MDAnalysisTests/data/lammps/albite_triclinic.dump
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
ITEM: TIMESTEP
hmacdope marked this conversation as resolved.
Show resolved Hide resolved
0
ITEM: NUMBER OF ATOMS
17
ITEM: BOX BOUNDS xy xz yz pp pp pp
-6.5875693349397721e+00 1.8337813315377392e+01 1.5067439154787670e+00
-5.4551678251502922e-01 2.5958964273996141e+01 -6.2664145519294436e+00
-4.5447071698045266e-02 1.2993982724334792e+01 -4.2179319547892025e-01
ITEM: ATOMS id type xs ys zs
192 1 0.204242 0.016215 0.0425371
85 1 0.154896 0.128649 0.358162
295 1 0.068964 0.161397 0.30523
300 1 0.231516 0.137347 0.326739
188 1 0.380162 0.0231137 0.125236
191 1 0.284665 0.224339 0.143678
299 1 0.122607 0.0686822 0.320504
159 1 0.16712 0.0505129 0.186929
136 1 0.476733 0.105929 0.221843
146 1 0.465467 0.020446 0.351053
193 1 0.212895 0.145535 0.00193565
81 1 0.0461335 0.227371 0.304292
189 1 0.375314 0.159839 0.0576947
43 1 0.438163 0.0223986 0.0606575
304 1 0.419694 0.108666 0.424637
86 1 0.163581 0.0158867 0.401886
302 1 0.363831 0.212981 0.38212
4 changes: 4 additions & 0 deletions testsuite/MDAnalysisTests/datafiles.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,10 +136,12 @@
"LAMMPScnt", "LAMMPScnt2", # triclinic box
"LAMMPShyd", "LAMMPShyd2",
"LAMMPSdata_deletedatoms", # with deleted atoms
"LAMMPSdata_triclinic", # lammpsdata file to test triclinic dimension parsing, albite with most atoms deleted
"LAMMPSDUMP",
"LAMMPSDUMP_long", # lammpsdump file with a few zeros sprinkled in the first column first frame
"LAMMPSDUMP_allcoords", # lammpsdump file with all coordinate conventions (x,xs,xu,xsu) present, from LAMMPS rdf example
"LAMMPSDUMP_nocoords", # lammpsdump file with no coordinates
"LAMMPSDUMP_triclinic", # lammpsdump file to test triclinic dimension parsing, albite with most atoms deleted
"unordered_res", # pdb file with resids non sequential
"GMS_ASYMOPT", # GAMESS C1 optimization
"GMS_SYMOPT", # GAMESS D4h optimization
Expand Down Expand Up @@ -490,10 +492,12 @@
LAMMPShyd = resource_filename(__name__, "data/lammps/hydrogen-class1.data")
LAMMPShyd2 = resource_filename(__name__, "data/lammps/hydrogen-class1.data2")
LAMMPSdata_deletedatoms = resource_filename(__name__, 'data/lammps/deletedatoms.data')
LAMMPSdata_triclinic = resource_filename(__name__, "data/lammps/albite_triclinic.data")
LAMMPSDUMP = resource_filename(__name__, "data/lammps/wat.lammpstrj.bz2")
LAMMPSDUMP_long = resource_filename(__name__, "data/lammps/wat.lammpstrj_long.bz2")
LAMMPSDUMP_allcoords = resource_filename(__name__, "data/lammps/spce_all_coords.lammpstrj.bz2")
LAMMPSDUMP_nocoords = resource_filename(__name__, "data/lammps/spce_no_coords.lammpstrj.bz2")
LAMMPSDUMP_triclinic = resource_filename(__name__, "data/lammps/albite_triclinic.dump")


unordered_res = resource_filename(__name__, "data/unordered_res.pdb")
Expand Down