Adding an XTC reporter (#4001)

* Preliminary work on XTC reporter
1. Move and adapt xtc writer/reader from moleculekit (explicit permission
granted by the authors to do so)
2. Create XTCTrajectoryFile
3. Create XTCReporter

* Add licence and attribution to c++ xtc library
Apply clang-format to it
Remove some unused functions and document the rest

* Add attribution and licence to cython wrappers for the xtc library
Remove some unused functions

* Change XTCTrajectoryFile to XTCFile
Simplify the interface and document the class

* Add test for the xtc file parser

* Update XTC reporter with new parser name

* Fix incorrect function name in XTCReporter

* XTCFile:
	* Add function to get number of frames
	* Add function to read a group of frames from a file
	* Add tests for the above

* Ensure data is passed as float32 in XTC file

* Add XTCReporter and tests

* Add more tests to XTCReporter

* Remove unnecessary pdb reporter in XTC tests

* Copy test xtc file in python/tests/systems to build directory for testing

* Remove XTC file reading from the interface
Make XTCFile mimic DCDFile more closely

* Use xtc_read to test the correctness of the XTC reporter

* Add a test for reporting triclinic boxes

* Make XTC library compatible with triclinic boxes.
Adapt XTCFile to triclinic boxes

* Change XTCFile to take a file as argument instead of a filename

* Match DCDFile handling of the box

* Fix comment

* Revert "Change XTCFile to take a file as argument instead of a filename"

This reverts commit 9815d4790b3886cc8a741586792268e80a227ba0.

* Fix dangling file name issue

* Remove index file functionality from XTC parser.
Remove unused define switch PLATFORM_Linux

* Fix formatting

* Remove inconsistent variable naming in xtcfile.py

* Change file argument name to match other reporters

* Do not turn off error checking in cython wrappers

* Fix leftover fileName in reporter

* Rewrite wrapper to xtclib in C++

* Small changes to wrapper code

* Small changes to wrapper code

* Small changes to wrapper code

* XTCFile: Get number of atoms directly from topology

* DCDFile: Get number of atoms directly from topology

* Change constexpr to const

* Check precision in XTC file matches the written one

* Add a write function to XTCFrame.
Make write check for errors C++ side.

* Rewrite large trajectory files without loading the whole file to memory

* Remove unused code in XTC test

* Avoid spurious copy of the positions array when calling xtc_write_frame

* Pass box as reference

* Remove unnecessary imports and definitions

* Fix formatting

* Use std::string instead of char*

* Use .c_str()  instead of .data()

* Fix crash in Mac by correctly checking precision

* Use TemporaryDirectory for tests instead of NamedTemporaryFile (Fixes windows ci)

* Remove unnecessary file creation

* Propagate exceptions via cython

* Switch to TemporaryDirectory in xtcfile.py

* Remove unnecessary include

* Update some comments and document functions

* Add XTC reporter to the docs

docs-source/usersguide/application/02_running_sims.rst

@@ -1379,13 +1379,14 @@ Writing Trajectories
 ====================
 
 
-OpenMM can save simulation trajectories to disk in three formats: PDB_,
-`PDBx/mmCIF`_, and DCD_.  All of these are widely supported formats, so you
+OpenMM can save simulation trajectories to disk in four formats: PDB_,
+`PDBx/mmCIF`_, DCD_ and XTC_.  All of these are widely supported formats, so you
 should be able to read them into most analysis and visualization programs.
 
 .. _PDB: http://www.wwpdb.org/documentation/format33/v3.3.html
 .. _PDBx/mmCIF: http://mmcif.wwpdb.org
 .. _DCD: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html
+.. _XTC: https://manual.gromacs.org/archive/5.0.4/online/xtc.html
 
 To save a trajectory, just add a “reporter” to the simulation, as shown in the
 example scripts above:
@@ -1394,9 +1395,9 @@ example scripts above:
     simulation.reporters.append(PDBReporter('output.pdb', 1000))
 
 The two parameters of the :class:`PDBReporter` are the output filename and how often (in
-number of time steps) output structures should be written.  To use PDBx/mmCIF or
-DCD format, just replace :class:`PDBReporter` with :class:`PDBxReporter` or
-:class:`DCDReporter`.  The parameters represent the same values:
+number of time steps) output structures should be written.  To use PDBx/mmCIF,
+DCD or XTC format, just replace :class:`PDBReporter` with :class:`PDBxReporter`, 
+:class:`DCDReporter` or :class:`XTCReporter`.  The parameters represent the same values:
 ::
 
     simulation.reporters.append(DCDReporter('output.dcd', 1000))

wrappers/python/CMakeLists.txt

@@ -51,10 +51,14 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*README.txt"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.py"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pyx"
+	"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pxd"
+	"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.cpp"
+	"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.h"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.i"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ini"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xtc"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdbx"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.cif"

wrappers/python/openmm/app/__init__.py

@@ -14,10 +14,12 @@ __email__ = "peastman@stanford.edu"
 
 from .topology import Topology, Chain, Residue, Atom
 from .pdbfile import PDBFile
+from .xtcfile import XTCFile
 from .pdbxfile import PDBxFile
 from .forcefield import ForceField
 from .simulation import Simulation
 from .pdbreporter import PDBReporter, PDBxReporter
+from .xtcreporter import XTCReporter
 from .amberprmtopfile import AmberPrmtopFile, HCT, OBC1, OBC2, GBn, GBn2
 from .amberinpcrdfile import AmberInpcrdFile
 from .dcdfile import DCDFile
@@ -53,4 +55,3 @@ Double = topology.Double
 Triple = topology.Triple
 Aromatic = topology.Aromatic
 Amide = topology.Amide
-

wrappers/python/openmm/app/dcdfile.py

@@ -117,7 +117,7 @@ class DCDFile(object):
         periodicBoxVectors : tuple of Vec3=None
             The vectors defining the periodic box.
         """
-        if len(list(self._topology.atoms())) != len(positions):
+        if self._topology.getNumAtoms() != len(positions):
             raise ValueError('The number of positions must match the number of atoms')
         if is_quantity(positions):
             positions = positions.value_in_unit(nanometers)

wrappers/python/openmm/app/internal/xtc_utils/include/xdrfile_xtc.h

+/* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- 
+ *
+ * $Id: xdrfile_xtc.h,v 1.5 2009/05/18 09:06:38 spoel Exp $
+ *
+ * Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
+ * Coordinate compression (c) by Frans van Hoesel. 
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 3
+ * of the License, or (at your option) any later version.
+ */
+
+#ifndef _xdrfile_xtc_h
+#define _xdrfile_xtc_h
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "xdrfile.h"
+  
+  /* All functions return exdrOK if succesfull. 
+   * (error codes defined in xdrfile.h).
+   */  
+   
+  /* This function returns the number of atoms in the xtc file in *natoms */
+  extern int read_xtc_natoms(char *fn,int *natoms);
+  
+  /* Read one frame of an open xtc file */
+  extern int read_xtc(XDRFILE *xd,int natoms,int *step,float *time,
+		      matrix box,rvec *x,float *prec);
+  
+  /* Write a frame to xtc file */
+  extern int write_xtc(XDRFILE *xd,
+		       int natoms,int step,float time,
+		       matrix box,rvec *x,float prec);
+  
+#ifdef __cplusplus
+}
+#endif
+
+#endif

wrappers/python/openmm/app/internal/xtc_utils/include/xtc.h

+/*
+MIT License
+
+Copyright (c) 2023 Accellera
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+Contributors: Stefan Doerr, Raul P. Pelaez
+*/
+#ifndef XTC
+#define XTC
+#include "xdrfile.h"
+#include<string>
+// Get the number of frames in a trajectory file
+int xtc_nframes(std::string filename);
+
+// Get the number of atoms in a trajectory file
+int xtc_natoms(std::string filename);
+
+// Reads nframes from  an xtc file, fills the arrays  in the input with
+// the data. Each array must be allocated with the correct size.
+void xtc_read(std::string filename, float* coords_arr, float* box_arr, float* time_arr, int* step_arr, int natoms, int nframes);
+
+// Appends a trajectory to a file.
+void xtc_write(std::string filename, int natoms, int nframes, int* step, float* timex, float* pos, float* box);
+
+// Rewrites a trajectory file with a new timestep and starting step number.
+// Useful when the step number is larger than 2^32.
+void xtc_rewrite_with_new_timestep(std::string filename_in, std::string filename_out, int first_step, int interval, float dt);
+#endif

wrappers/python/openmm/app/internal/xtc_utils/src/xdrfile_xtc.cpp

+/* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
+ *
+ * $Id: xdrfile_xtc.c,v 1.5 2009/05/18 09:06:38 spoel Exp $
+ *
+ * Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
+ * Coordinate compression (c) by Frans van Hoesel.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 3
+ * of the License, or (at your option) any later version.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include "xdrfile.h"
+#include "xdrfile_xtc.h"
+
+#define MAGIC 1995
+enum { FALSE, TRUE };
+
+static int xtc_header(XDRFILE *xd,int *natoms,int *step,float *time,mybool bRead)
+{
+	int result,magic,n=1;
+
+	/* Note: read is same as write. He he he */
+	magic  = MAGIC;
+	if ((result = xdrfile_write_int(&magic,n,xd)) != n)
+		{
+			if (bRead)
+				return exdrENDOFFILE;
+			else
+				return exdrINT;
+		}
+	if (magic != MAGIC)
+		return exdrMAGIC;
+	if ((result = xdrfile_write_int(natoms,n,xd)) != n)
+		return exdrINT;
+	if ((result = xdrfile_write_int(step,n,xd)) != n)
+		return exdrINT;
+	if ((result = xdrfile_write_float(time,n,xd)) != n)
+		return exdrFLOAT;
+
+	return exdrOK;
+}
+
+static int xtc_coord(XDRFILE *xd,int *natoms,matrix box,rvec *x,float *prec,
+					 mybool bRead)
+{
+	int result;
+
+	/* box */
+	result = xdrfile_read_float(box[0],DIM*DIM,xd);
+	if (DIM*DIM != result)
+		return exdrFLOAT;
+	else
+		{
+			if (bRead)
+				{
+					result = xdrfile_decompress_coord_float(x[0],natoms,prec,xd);
+					if (result != *natoms)
+						return exdr3DX;
+				}
+			else
+				{
+					result = xdrfile_compress_coord_float(x[0],*natoms,*prec,xd);
+					if (result != *natoms)
+						return exdr3DX;
+				}
+		}
+	return exdrOK;
+}
+
+int read_xtc_natoms(char *fn,int *natoms)
+{
+	XDRFILE *xd;
+	int step,result;
+	float time;
+
+	xd = xdrfile_open(fn,"r");
+	if (NULL == xd)
+		return exdrFILENOTFOUND;
+	result = xtc_header(xd,natoms,&step,&time,TRUE);
+	xdrfile_close(xd);
+
+	return result;
+}
+
+int read_xtc(XDRFILE *xd,
+			 int natoms,int *step,float *time,
+			 matrix box,rvec *x,float *prec)
+/* Read subsequent frames */
+{
+	int result;
+
+	if ((result = xtc_header(xd,&natoms,step,time,TRUE)) != exdrOK)
+		return result;
+
+	if ((result = xtc_coord(xd,&natoms,box,x,prec,1)) != exdrOK)
+		return result;
+
+	return exdrOK;
+}
+
+int write_xtc(XDRFILE *xd,
+			  int natoms,int step,float time,
+			  matrix box,rvec *x,float prec)
+/* Write a frame to xtc file */
+{
+	int result;
+
+	if ((result = xtc_header(xd,&natoms,&step,&time,FALSE)) != exdrOK)
+		return result;
+
+	if ((result = xtc_coord(xd,&natoms,box,x,&prec,0)) != exdrOK)
+		return result;
+
+	return exdrOK;
+}

wrappers/python/openmm/app/internal/xtc_utils/src/xtc.cpp

+/*
+MIT License
+
+Copyright (c) 2023 Accellera
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+Contributors: Stefan Doerr, Raul P. Pelaez
+*/
+#include "xtc.h"
+#include "xdrfile_xtc.h"
+#include <vector>
+#include <string>
+#include <stdexcept>
+
+// Helper functions to convert between atom/frame and x/y/z indices
+static size_t Xf(size_t atom, size_t frame, size_t nframes) {
+    return atom * 3 * nframes + frame;
+}
+
+static size_t Yf(size_t xidx, size_t nframes) {
+    return xidx + nframes;
+}
+
+static size_t Zf(size_t yidx, size_t nframes) {
+    return yidx + nframes;
+}
+
+// Helper struct to manage the XDRFILE pointer
+struct XDRFILE_RAII {
+    XDRFILE* xd;
+    XDRFILE_RAII(std::string filename, std::string mode) : xd(xdrfile_open(filename.c_str(), mode.c_str())) {
+        if (!xd) {
+            throw std::runtime_error("xtc file: Could not open file");
+        }
+    }
+
+    ~XDRFILE_RAII() {
+        if (xd)
+            xdrfile_close(xd);
+    }
+
+    operator XDRFILE*() const {
+        return xd;
+    }
+};
+
+int xtc_natoms(std::string filename) {
+    int natoms = 0;
+    if (exdrOK != read_xtc_natoms(const_cast<char*>(filename.c_str()), &natoms)) {
+        throw std::runtime_error("xtc_read(): could not get natoms\n");
+    }
+    return natoms;
+}
+
+struct XTCFrame {
+    int step;
+    float time;
+    matrix box;
+    std::vector<float> positions;
+    int natoms;
+    const float prec = 1000.0;
+    XTCFrame(int natoms) : positions(3 * natoms), natoms(natoms) {
+    }
+
+    // Read the next frame from the XTC file and store it in this object
+    int readNextFrame(XDRFILE* xd) {
+        float in_prec;
+        auto* p_ptr = reinterpret_cast<rvec*>(positions.data());
+        int status = read_xtc(xd, natoms, &step, &time, box, p_ptr, &in_prec);
+        if (status == exdrOK && prec != in_prec) {
+            throw std::runtime_error("xtc_read(): precision mismatch\n");
+        }
+        if (status == exdr3DX) {
+            throw std::runtime_error("xtc_read(): XTC file is corrupt\n");
+        }
+        return status;
+    }
+
+    // Write the current frame to the XTC file
+    int appendFrameToFile(XDRFILE* xd) {
+        auto* p_ptr = reinterpret_cast<rvec*>(positions.data());
+        int err = write_xtc(xd, natoms, step, time, box, p_ptr, prec);
+        if (err != exdrOK) {
+            throw std::runtime_error("xtc_write(): could not write frame\n");
+        }
+        return err;
+    }
+};
+
+int xtc_nframes(std::string filename) {
+    int nframes = 0;
+    int natoms = xtc_natoms(filename);
+    if (!natoms) {
+        throw std::runtime_error("xtc_read(): natoms is 0\n");
+    }
+    XDRFILE_RAII xd(filename, "r");
+    XTCFrame frame(natoms);
+    while (exdrOK == frame.readNextFrame(xd)) {
+        nframes++;
+    }
+    return nframes;
+}
+
+void xtc_read(std::string filename, float* coords_arr, float* box_arr, float* time_arr, int* step_arr, int natoms, int nframes) {
+    if (natoms == 0) {
+        throw std::runtime_error("xtc_read(): natoms is 0\n");
+    }
+    XDRFILE_RAII xd(filename, "r");
+    int fidx = 0;
+    XTCFrame frame(natoms);
+    while (exdrOK == frame.readNextFrame(xd)) {
+        time_arr[fidx] = frame.time;
+        step_arr[fidx] = frame.step;
+        for (int i = 0; i < 3; i++) {
+            for (int j = 0; j < 3; j++) {
+                box_arr[fidx + (3 * i + j) * nframes] = frame.box[i][j];
+            }
+        }
+        for (int aidx = 0; aidx < natoms; aidx++) {
+            int xidx = Xf(aidx, fidx, nframes);
+            int yidx = Yf(xidx, nframes);
+            int zidx = Zf(yidx, nframes);
+            coords_arr[xidx] = frame.positions[3 * aidx + 0];
+            coords_arr[yidx] = frame.positions[3 * aidx + 1];
+            coords_arr[zidx] = frame.positions[3 * aidx + 2];
+        }
+        fidx++;
+    }
+}
+
+static void box_from_array(matrix& matrix_box, float* box, int frame, int nframes) {
+    for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+            matrix_box[i][j] = box[(3 * i + j) * nframes + frame];
+        }
+    }
+}
+
+void xtc_write(std::string filename, int natoms, int nframes, int* step, float* timex, float* pos, float* box) {
+    XDRFILE_RAII xd(filename, "a");
+    XTCFrame frame(natoms);
+    for (int f = 0; f < nframes; f++) {
+        box_from_array(frame.box, box, f, nframes);
+        for (int i = 0; i < natoms; i++) {
+            int xidx = Xf(i, f, nframes);
+            int yidx = Yf(xidx, nframes);
+            int zidx = Zf(yidx, nframes);
+            frame.positions[3 * i + 0] = pos[xidx];
+            frame.positions[3 * i + 1] = pos[yidx];
+            frame.positions[3 * i + 2] = pos[zidx];
+        }
+        frame.step = step[f];
+        frame.time = timex[f];
+        frame.appendFrameToFile(xd);
+    }
+}
+
+void xtc_rewrite_with_new_timestep(std::string filename_in, std::string filename_out, int first_step, int interval, float dt) {
+    int natoms = xtc_natoms(filename_in);
+    if (natoms == 0) {
+        throw std::runtime_error("xtc_read(): natoms is 0\n");
+    }
+    XDRFILE_RAII xd_in(filename_in, "r");
+    XDRFILE_RAII xd_out(filename_out, "a");
+    XTCFrame frame(natoms);
+    int i = 0;
+    while (exdrOK == frame.readNextFrame(xd_in)) {
+        frame.step = first_step + i * interval;
+        frame.time = frame.step * dt;
+        frame.appendFrameToFile(xd_out);
+        i++;
+    }
+}

wrappers/python/openmm/app/internal/xtc_utils/xtc.pyx

+# MIT License
+
+# Copyright (c) 2023 Accellera
+
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+# Contributors: Stefan Doerr, Raul P. Pelaez
+
+import numpy as np
+cimport numpy as np
+cimport xtclib
+from libcpp.string cimport string
+ctypedef np.float32_t FLOAT32_t
+
+def get_xtc_nframes(string filename):
+    """
+    Get the number of frames in a xtc file.
+    Parameters
+    ----------
+    filename: string
+        The filename of the xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    Returns
+    -------
+    nframes: int
+        The number of frames in the xtc file
+    """
+    return xtclib.xtc_nframes(filename)
+
+def get_xtc_natoms(string filename):
+    """
+    Get the number of atoms in a xtc file.
+    Parameters
+    ----------
+    filename: string
+        The filename of the xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    Returns
+    -------
+    natoms: int
+        The number of atoms in the xtc file
+    """
+    return xtclib.xtc_natoms(filename)
+
+def read_xtc(string filename):
+    """
+    Reads a xtc file and return its contents.
+
+    Parameters
+    ----------
+    filename: string
+        The filename of the xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    Returns
+    -------
+    coords: np.ndarray
+        The coordinates of the atoms in the xtc file. Shape: (n_atoms, 3, n_frames)
+    box: np.ndarray
+        The box vectors of the xtc file. Shape: (3, 3, n_frames)
+    time: np.ndarray
+        The time of each frame. Shape: (n_frames,)
+    step: np.ndarray
+        The step of each frame. Shape: (n_frames,)
+    """
+    cdef int natoms = get_xtc_natoms(filename)
+    cdef int nframes = get_xtc_nframes(filename)
+
+    cdef FLOAT32_t[:, :, ::1] coords = np.zeros((natoms, 3, nframes), dtype=np.float32)
+    cdef FLOAT32_t[:, :, ::1] box = np.zeros((3, 3, nframes), dtype=np.float32)
+    cdef FLOAT32_t[::1] time = np.zeros(nframes, dtype=np.float32)
+    cdef int[::1] step = np.zeros(nframes, dtype=np.int32)
+
+    xtclib.xtc_read(
+        filename,
+        &coords[0, 0, 0],
+        &box[0, 0, 0],
+        &time[0],
+        &step[0],
+        natoms,
+        nframes,
+    )
+    return np.asarray(coords), np.asarray(box), np.asarray(time), np.asarray(step)
+
+def xtc_write_frame(string filename, float[:, :] coords, float[:, :] box, float time, int step):
+    """
+    Appends a single frame to a xtc file (if the file does not exist it is created by this function).
+    Parameters
+    ----------
+    filename: string
+        The filename of the xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    coords: np.ndarray
+        The coordinates of the atoms in the frame. Shape: (n_atoms, 3)
+    box: np.ndarray
+        The box vectors of the frame. Shape: (3, 3)
+    time: float
+        The time of the frame
+    step: int
+        The step of the frame
+    """
+    cdef int natoms = coords.shape[0]
+    cdef int nframes = 1
+    xtclib.xtc_write(
+        filename,
+        natoms,
+        nframes,
+        &step,
+        &time,
+        &coords[0, 0],
+        &box[0, 0]
+    )
+
+
+def xtc_rewrite_with_new_timestep(string filename_in, string filename_out,
+				  int first_step, int interval, float dt):
+    """
+    Rewrites a trajectory file with a new timestep and starting step number.
+    Parameters
+    ----------
+    filename_in: string
+        The filename of the input xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    filename_out: string
+        The filename of the output xtc file. You need to pass the string with filename.encode("UTF-8") to this function
+    first_step: int
+        The first step to be written to the output file
+    interval: int
+        The interval between steps to be written to the output file
+    dt: float
+        The timestep of the output file
+    """
+    xtclib.xtc_rewrite_with_new_timestep(filename_in, filename_out, first_step, interval, dt)

wrappers/python/openmm/app/internal/xtc_utils/xtclib.pxd

+# MIT License
+
+# Copyright (c) 2023 Accellera
+
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+# Contributors: Stefan Doerr, Raul P. Pelaez
+from libcpp.string cimport string
+cdef extern from "include/xtc.h":
+    cdef int xtc_nframes(string filename) except +
+    cdef int xtc_natoms(string filename) except +
+    cdef void xtc_read(string filename, float *coords_arr, float *box_arr, float *time_arr, int *step_arr, int natoms, int nframes) except +
+    cdef void xtc_write(string filename, int natoms, int nframes, int *step, float *timex, float *pos, float *box) except +
+    cdef void xtc_rewrite_with_new_timestep(string filename_in, string filename_out,
+				  int first_step, int interval, float dt) except +

wrappers/python/openmm/app/xtcfile.py

+from __future__ import print_function, division, absolute_import
+
+__author__ = "Raul P. Pelaez"
+__version__ = "1.0"
+
+from openmm.app.internal.xtc_utils import (
+    xtc_rewrite_with_new_timestep,
+    xtc_write_frame,
+    get_xtc_nframes,
+    get_xtc_natoms,
+)
+import numpy as np
+import os
+from openmm import Vec3
+from openmm.unit import nanometers, femtoseconds, is_quantity, norm
+import math
+import tempfile
+import shutil
+
+
+class XTCFile(object):
+
+    """XTCFile provides methods for creating XTC files.
+    To use this class, create a XTCFile object, then call writeModel() once for each model in the file.
+    """
+
+    def __init__(self, fileName, topology, dt, firstStep=0, interval=1, append=False):
+        """Create a XTC file, or open an existing file to append.
+
+        Parameters
+        ----------
+        fileName : str
+            A file name to write to
+        topology : Topology
+            The Topology defining the molecular system being written
+        dt : time
+            The time step used in the trajectory
+        firstStep : int=0
+            The index of the first step in the trajectory
+        interval : int=1
+            The frequency (measured in time steps) at which states are written
+            to the trajectory
+        append : bool=False
+            If True, open an existing XTC file to append to.  If False, create a new file.
+        """
+        if not isinstance(fileName, str):
+            raise TypeError("fileName must be a string")
+        self._filename = fileName
+        self._topology = topology
+        self._firstStep = firstStep
+        self._interval = interval
+        self._modelCount = 0
+        if is_quantity(dt):
+            dt = dt.value_in_unit(femtoseconds)
+        self._dt = dt
+        if append:
+            if not os.path.isfile(self._filename):
+                raise FileNotFoundError(f"The file '{self._filename}' does not exist.")
+            self._modelCount = get_xtc_nframes(self._filename.encode("utf-8"))
+            natoms = get_xtc_natoms(self._filename.encode("utf-8"))
+            if natoms != len(list(topology.atoms())):
+                raise ValueError(
+                    f"The number of atoms in the topology ({len(list(topology.atoms()))}) does not match the number of atoms in the XTC file ({natoms})"
+                )
+        else:
+            if os.path.isfile(self._filename) and os.path.getsize(self._filename) > 0:
+                raise FileExistsError(f"The file '{self._filename}' already exists.")
+
+    def _getNumFrames(self):
+        return get_xtc_nframes(self._filename.encode("utf-8"))
+
+    def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None):
+        """Write out a model to the XTC file.
+
+        The periodic box can be specified either by the unit cell dimensions
+        (for a rectangular box), or the full set of box vectors (for an
+        arbitrary triclinic box).  If neither is specified, the box vectors
+        specified in the Topology will be used. Regardless of the value
+        specified, no dimensions will be written if the Topology does not
+        represent a periodic system.
+
+        Parameters
+        ----------
+        positions : list
+            The list of atomic positions to write
+        unitCellDimensions : Vec3=None
+            The dimensions of the crystallographic unit cell.
+        periodicBoxVectors : tuple of Vec3=None
+            The vectors defining the periodic box.
+        """
+        if self._topology.getNumAtoms() != len(positions):
+            raise ValueError("The number of positions must match the number of atoms")
+        if is_quantity(positions):
+            positions = positions.value_in_unit(nanometers)
+        if any(math.isnan(norm(pos)) for pos in positions):
+            raise ValueError(
+                "Particle position is NaN.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan"
+            )
+        if any(math.isinf(norm(pos)) for pos in positions):
+            raise ValueError(
+                "Particle position is infinite.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan"
+            )
+
+        self._modelCount += 1
+        if (
+            self._interval > 1
+            and self._firstStep + self._modelCount * self._interval > 1 << 31
+        ):
+            # This will exceed the range of a 32 bit integer.  To avoid crashing or producing a corrupt file,
+            # update the file to say the trajectory consisted of a smaller number of larger steps (so the
+            # total trajectory length remains correct).
+            self._firstStep //= self._interval
+            self._dt *= self._interval
+            self._interval = 1
+            with tempfile.TemporaryDirectory() as temp:
+                fname = os.path.join(temp, "temp.xtc")
+                xtc_rewrite_with_new_timestep(
+                    self._filename.encode("utf-8"),
+                    fname.encode("utf-8"),
+                    self._firstStep,
+                    self._interval,
+                    self._dt,
+                )
+                shutil.copyfile(fname, self._filename)
+        boxVectors = self._topology.getPeriodicBoxVectors()
+        if boxVectors is not None:
+            if periodicBoxVectors is not None:
+                boxVectors = periodicBoxVectors
+                if is_quantity(boxVectors):
+                    boxVectors = boxVectors.value_in_unit(nanometers)
+            elif unitCellDimensions is not None:
+                if is_quantity(unitCellDimensions):
+                    unitCellDimensions = unitCellDimensions.value_in_unit(nanometers)
+                boxVectors = (
+                    Vec3(unitCellDimensions[0], 0, 0),
+                    Vec3(0, unitCellDimensions[1], 0),
+                    Vec3(0, 0, unitCellDimensions[2]),
+                )
+            boxVectors = np.array(
+                [[vec.x, vec.y, vec.z] for vec in boxVectors], dtype=np.float32
+            )
+        else:
+            boxVectors = np.zeros((3, 3)).astype(np.float32)
+        step = self._modelCount * self._interval + self._firstStep
+        time = step * self._dt
+        xtc_write_frame(
+            self._filename.encode("utf-8"),
+            np.array(positions, dtype=np.float32),
+            boxVectors,
+            np.float32(time),
+            np.int32(step),
+        )

wrappers/python/openmm/app/xtcreporter.py

+__author__ = "Raul P. Pelaez"
+from openmm.app import XTCFile
+
+class XTCReporter(object):
+    """XTCReporter outputs a series of frames from a Simulation to a XTC file.
+
+    To use it, create a XTCReporter, then add it to the Simulation's list of reporters.
+    """
+
+    def __init__(self, file, reportInterval, append=False, enforcePeriodicBox=None):
+        """Create a XTCReporter.
+
+        Parameters
+        ----------
+        file : string
+            The file to write to
+        reportInterval : int
+            The interval (in time steps) at which to write frames
+        append : bool=False
+            If True, open an existing XTC file to append to.  If False, create a new file.
+        enforcePeriodicBox: bool
+            Specifies whether particle positions should be translated so the center of every molecule
+            lies in the same periodic box.  If None (the default), it will automatically decide whether
+            to translate molecules based on whether the system being simulated uses periodic boundary
+            conditions.
+        """
+        self._reportInterval = reportInterval
+        self._append = append
+        self._enforcePeriodicBox = enforcePeriodicBox
+        self._fileName = file
+        self._xtc = None
+
+    def describeNextReport(self, simulation):
+        """Get information about the next report this object will generate.
+
+        Parameters
+        ----------
+        simulation : Simulation
+            The Simulation to generate a report for
+
+        Returns
+        -------
+        tuple
+            A six element tuple. The first element is the number of steps
+            until the next report. The next four elements specify whether
+            that report will require positions, velocities, forces, and
+            energies respectively.  The final element specifies whether
+            positions should be wrapped to lie in a single periodic box.
+        """
+        steps = self._reportInterval - simulation.currentStep % self._reportInterval
+        return (steps, True, False, False, False, self._enforcePeriodicBox)
+
+    def report(self, simulation, state):
+        """Generate a report.
+
+        Parameters
+        ----------
+        simulation : Simulation
+            The Simulation to generate a report for
+        state : State
+            The current state of the simulation
+        """
+
+        if self._xtc is None:
+            self._xtc = XTCFile(
+                self._fileName,
+                simulation.topology,
+                simulation.integrator.getStepSize(),
+                simulation.currentStep,
+                self._reportInterval,
+                self._append,
+            )
+        self._xtc.writeModel(
+            state.getPositions(), periodicBoxVectors=state.getPeriodicBoxVectors()
+        )

wrappers/python/setup.py

@@ -229,6 +229,22 @@ def buildKeywordDictionary(major_version_num=MAJOR_VERSION_NUM,
     setupKeywords["ext_modules"] = [Extension(**extensionArgs)]
     setupKeywords["ext_modules"] += cythonize('openmm/app/internal/*.pyx')
 
+    setupKeywords["ext_modules"] +=cythonize(Extension(
+        "openmm.app.internal.xtc_utils",
+        sources=[
+            "openmm/app/internal/xtc_utils/src/xdrfile_xtc.cpp",
+            "openmm/app/internal/xtc_utils/src/xdrfile.cpp",
+            "openmm/app/internal/xtc_utils/src/xtc.cpp",
+            "openmm/app/internal/xtc_utils/xtc.pyx",
+        ],
+        include_dirs=include_dirs +[
+            "openmm/app/internal/xtc_utils/include",
+            "openmm/app/internal/xtc_utils/",
+            numpy.get_include(),
+        ],
+        language="c++",
+    ))
+
     outputString = ''
     firstTab     = 40
     secondTab    = 60

wrappers/python/tests/TestXtcFile.py

+__author__ = "Raul P. Pelaez"
+import sys
+import os
+import unittest
+from openmm import unit
+import tempfile
+from openmm import app
+from random import random
+import openmm as mm
+import numpy as np
+from openmm.app.internal.xtc_utils import read_xtc
+
+class TestXtcFile(unittest.TestCase):
+    def test_xtc_triclinic(self):
+        """Test the XTC file by writing a trajectory and reading it back. Using a triclinic box"""
+        with tempfile.TemporaryDirectory() as temp:
+            fname = os.path.join(temp, 'traj.xtc')
+            pdbfile = app.PDBFile("systems/alanine-dipeptide-implicit.pdb")
+            # Set some arbitrary size for the unit cell so that a box is included in the trajectory
+            pdbfile.topology.setUnitCellDimensions([10, 10, 10])
+            natom = len(list(pdbfile.topology.atoms()))
+            nframes = 20
+            xtc = app.XTCFile(fname, pdbfile.topology, 0.001)
+            coords = []
+            box = []
+            for i in range(nframes):
+                coords.append(
+                    [mm.Vec3(random(), random(), random()) for j in range(natom)]
+                    * unit.nanometers
+                )
+                box_i = (
+                    mm.Vec3(random(), random(), random()) * unit.nanometers,
+                    mm.Vec3(random(), random(), random()) * unit.nanometers,
+                    mm.Vec3(random(), random(), random()) * unit.nanometers,
+                )
+                box.append(np.array([[vec.x, vec.y, vec.z] for vec in box_i]))
+                xtc.writeModel(coords[i], periodicBoxVectors=box_i)
+            # The  XTCFile class  does not  provide a  way to  read the
+            # trajectory back, but the underlying XTC library does
+            coords_read, box_read, time, step = read_xtc(fname.encode("utf-8"))
+            self.assertEqual(coords_read.shape, (natom, 3, nframes))
+            self.assertEqual(box_read.shape, (3, 3, nframes))
+            self.assertEqual(len(time), nframes)
+            self.assertEqual(len(step), nframes)
+            coords = np.array(
+                [c.value_in_unit(unit.nanometers) for c in coords]
+            ).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(coords_read, coords, atol=1e-3))
+            box = np.array(box).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(box_read, box, atol=1e-3))
+            self.assertTrue(
+                np.allclose(time, np.arange(1, nframes + 1) * 0.001, atol=1e-5)
+            )
+            self.assertTrue(np.allclose(step, np.arange(1, nframes + 1), atol=1e-5))
+
+    def test_xtc_cubic(self):
+        """Test the XTC file by writing a trajectory and reading it back. Using a cubic box"""
+        with tempfile.TemporaryDirectory() as temp:
+            fname = os.path.join(temp, 'traj.xtc')
+            pdbfile = app.PDBFile("systems/alanine-dipeptide-implicit.pdb")
+            # Set some arbitrary size for the unit cell so that a box is included in the trajectory
+            pdbfile.topology.setUnitCellDimensions([10, 10, 10])
+            natom = len(list(pdbfile.topology.atoms()))
+            nframes=20
+            xtc = app.XTCFile(fname, pdbfile.topology, 0.001)
+            coords = []
+            box = []
+            for i in range(nframes):
+                coords.append(
+                    [mm.Vec3(random(), random(), random()) for j in range(natom)]
+                    * unit.nanometers
+                )
+                box_i = mm.Vec3(random(), random(), random()) * unit.nanometers
+                box.append(np.diag(box_i[:3].value_in_unit(unit.nanometers)))
+                xtc.writeModel(coords[i], unitCellDimensions=box_i)
+            #The  XTCFile class  does not  provide a  way to  read the
+            #trajectory back, but the underlying XTC library does
+            coords_read, box_read, time, step = read_xtc(fname.encode("utf-8"))
+            self.assertEqual(coords_read.shape, (natom,3,nframes))
+            self.assertEqual(box_read.shape, (3,3,nframes))
+            self.assertEqual(len(time), nframes)
+            self.assertEqual(len(step), nframes)
+            coords = np.array(
+                [c.value_in_unit(unit.nanometers) for c in coords]
+            ).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(coords_read, coords, atol=1e-3))
+            box = np.array(box).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(box_read, box, atol=1e-3))
+            self.assertTrue(
+                np.allclose(time, np.arange(1, nframes + 1) * 0.001, atol=1e-5)
+            )
+            self.assertTrue(np.allclose(step, np.arange(1, nframes + 1), atol=1e-5))
+
+    def test_xtc_box_from_topology(self):
+        """Test the XTC file by writing a trajectory and reading it back. Letting the box be set from the topology"""
+        with tempfile.TemporaryDirectory() as temp:
+            fname = os.path.join(temp, 'traj.xtc')
+            pdbfile = app.PDBFile("systems/alanine-dipeptide-implicit.pdb")
+            # Set some arbitrary size for the unit cell so that a box is included in the trajectory
+            unitCell = mm.Vec3(random(), random(), random()) * unit.nanometers
+            pdbfile.topology.setUnitCellDimensions(unitCell)
+            natom = len(list(pdbfile.topology.atoms()))
+            nframes = 20
+            xtc = app.XTCFile(fname, pdbfile.topology, 0.001)
+            coords = []
+            for i in range(nframes):
+                coords.append(
+                    [mm.Vec3(random(), random(), random()) for j in range(natom)]
+                    * unit.nanometers
+                )
+                xtc.writeModel(coords[i])
+            # The  XTCFile class  does not  provide a  way to  read the
+            # trajectory back, but the underlying XTC library does
+            coords_read, box_read, time, step = read_xtc(fname.encode("utf-8"))
+            self.assertEqual(coords_read.shape, (natom, 3, nframes))
+            self.assertEqual(box_read.shape, (3, 3, nframes))
+            self.assertEqual(len(time), nframes)
+            self.assertEqual(len(step), nframes)
+            coords = np.array(
+                [c.value_in_unit(unit.nanometers) for c in coords]
+            ).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(coords_read, coords, atol=1e-3))
+            boxVectors = (
+                mm.Vec3(unitCell[0].value_in_unit(unit.nanometers), 0, 0),
+                mm.Vec3(0, unitCell[1].value_in_unit(unit.nanometers), 0),
+                mm.Vec3(0, 0, unitCell[2].value_in_unit(unit.nanometers)),
+            )
+            boxVectors = np.array(
+                [[vec.x, vec.y, vec.z] for vec in boxVectors], dtype=np.float32
+            )
+            box = np.array(np.tile(boxVectors, (nframes, 1, 1))).transpose(1, 2, 0)
+            self.assertTrue(np.allclose(box_read, box, atol=1e-3))
+            self.assertTrue(
+                np.allclose(time, np.arange(1, nframes + 1) * 0.001, atol=1e-5)
+            )
+            self.assertTrue(np.allclose(step, np.arange(1, nframes + 1), atol=1e-5))
+
+    def testLongTrajectory(self):
+        """Test writing a trajectory that has more than 2^31 steps."""
+        with tempfile.TemporaryDirectory() as temp:
+            fname = os.path.join(temp, 'traj.xtc')
+            pdbfile = app.PDBFile("systems/alanine-dipeptide-implicit.pdb")
+            natom = len(list(pdbfile.topology.atoms()))
+            xtc = app.XTCFile(fname, pdbfile.topology, 0.001, interval=1000000000)
+            for i in range(5):
+                xtc.writeModel(
+                    [mm.Vec3(random(), random(), random()) for j in range(natom)]
+                    * unit.angstroms
+                )
+
+    def testAppend(self):
+        """Test appending to an existing trajectory."""
+        with tempfile.TemporaryDirectory() as temp:
+            fname = os.path.join(temp, 'traj.xtc')
+            pdb = app.PDBFile("systems/alanine-dipeptide-implicit.pdb")
+            ff = app.ForceField("amber99sb.xml", "tip3p.xml")
+            system = ff.createSystem(pdb.topology)
+
+            # Create a simulation and write some frames to a XTC file.
+
+            integrator = mm.VerletIntegrator(0.001 * unit.picoseconds)
+            simulation = app.Simulation(
+                pdb.topology,
+                system,
+                integrator,
+                mm.Platform.getPlatformByName("Reference"),
+            )
+            xtc = app.XTCReporter(fname, 2)
+            simulation.reporters.append(xtc)
+            simulation.context.setPositions(pdb.positions)
+            simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
+            simulation.step(10)
+            self.assertEqual(5, xtc._xtc._modelCount)
+            self.assertEqual(5, xtc._xtc._getNumFrames())
+            del simulation
+            del xtc
+
+            # Create a new simulation and have it append some more frames.
+
+            integrator = mm.VerletIntegrator(0.001 * unit.picoseconds)
+            simulation = app.Simulation(
+                pdb.topology,
+                system,
+                integrator,
+                mm.Platform.getPlatformByName("Reference"),
+            )
+            xtc = app.XTCReporter(fname, 2, append=True)
+            simulation.reporters.append(xtc)
+            simulation.context.setPositions(pdb.positions)
+            simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
+            simulation.step(10)
+            self.assertEqual(10, xtc._xtc._modelCount)
+            self.assertEqual(10, xtc._xtc._getNumFrames())
+            del simulation
+            del xtc
+
+
+if __name__ == "__main__":
+    unittest.main()