Implement Gromacs internal Whole Procedure

pygmx/__init__.py

@@ -7,7 +7,7 @@ Trajectories in trr format may be read with `.gromacs.reader.TRRReader`, which i
 
 import os
 
-from .tpxio import TPXReader
+from .tpxio import TPXReader, make_xtcframe_whole
 from .xtcio import XTCReader, read_xtcframe, append_xtcfile
 from .enxio import EDRFile
 from .errors import FileTypeError

pygmx/topology.pxd

@@ -68,6 +68,9 @@ cdef extern from "gromacs/topology/idef.h":
     ctypedef struct t_ilist:
         pass
 
+    ctypedef struct t_idef:
+        pass
+
     #cdef enum t_ft_enum:
     #    F_LJ
 
@@ -143,14 +146,45 @@ cdef extern from "gromacs/topology/topology.h":
         gmx_groups_t     groups
         t_symtab         symtab                      # The symbol table                     */
 
+    ctypedef struct t_topology:
+        char          **name                       # /* Name of the topology                 */
+        t_idef          idef                       # /* The interaction function definition  */
+        t_atoms         atoms                      # /* The atoms                            */
+        t_atomtypes     atomtypes                  # /* Atomtype properties                  */
+        t_block         cgs                        # /* The charge groups                    */
+        t_block         mols                       # /* The molecules                        */
+        gmx_bool        bIntermolecularInteractions# /* Inter.mol. int. ?   */
+        t_blocka        excls                      # /* The exclusions                       */
+        t_symtab        symtab                     # /* The symbol table                     */
+
+    ctypedef struct gmx_localtop_t:
+        t_idef        idef        # /* The interaction function definition  */
+        #t_atomtypes   atomtypes   # /* Atomtype properties                  */
+        #t_block       cgs         # /* The charge groups                    */
+        #   t_blocka      excls       # /* The exclusions                       */
+
+    void init_mtop(gmx_mtop_t *mtop)
+    void done_top(t_topology *top)
+
 # cdef extern from "gromacs/topology/topology.h":
         # generate a t_atoms struct for the system from gmx_mtop_t
         # t_atoms* mtop2atoms(gmx_mtop_t *mtop)
 cdef extern from "gromacs/topology/mtop_util.h":
     t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
+    t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bint freeMTop)
+    gmx_localtop_t *gmx_mtop_generate_local_top(const gmx_mtop_t *mtop, bint freeEnergyInteractionsAtEnd)
+
+
 
 cdef extern from "gromacs/pbcutil/rmpbc.h":
-    void rm_gropbc(const t_atoms *atoms, rvec x[], const matrix box);
+    ctypedef struct gmx_rmpbc_t:
+        pass
+
+    gmx_rmpbc_t gmx_rmpbc_init(const t_idef *idef, int ePBC, int natoms)
+    void gmx_rmpbc_done(gmx_rmpbc_t gpbc)
+    void rm_gropbc(const t_atoms *atoms, rvec x[], const matrix box)
+    void gmx_rmpbc(gmx_rmpbc_t gpbc, int natoms, const matrix box, rvec x[])
+
 
 
 

pygmx/tpxio.pyx

@@ -2,9 +2,11 @@
 # Cython wrapper around tpxio.cpp
 
 from libc cimport stdio
+from libc.string cimport memcpy
 
 import numpy as np
-#cimport numpy as np
+import cython
+cimport numpy as np
 
 from utility cimport *
 from math cimport *
@@ -254,9 +256,20 @@ cdef class TPXReader:
 
 
 @cython.binding(True)
-def make_xtcframe_whole(coords, box, <TPXReader>reader):
-    cdef t_atoms = gmx_mtop_global_atoms(reader.topology)
+def make_xtcframe_whole(coords, box, TPXReader reader):
+    
+    cdef int natoms = reader.topology.natoms
+    cdef gmx_localtop_t *top = gmx_mtop_generate_local_top(&reader.topology, True)
+    cdef gmx_rmpbc_t gpbc = gmx_rmpbc_init(&top.idef, -1, natoms)
+
     cdef np.ndarray[real, ndim=2] b = np.asarray(box, dtype=np.float32)
-    cdef np.ndarray[real, ndim=2] x = np.asarray(coords, dtype=np.float32)
-    rm_gropbc(const t_atoms *atoms, <rvec *>x.data, <matrix> b)
+    cdef np.ndarray[real, ndim=2] x = np.array(coords, dtype=np.float32).copy()
+    
+    gmx_rmpbc(gpbc, natoms, <rvec *>b.data, <rvec *>x.data)
+
+    # free up memory - in fact, we should free memory of 'top' too
+    # however, it shares memory with 'reader.topology', so we can not free it correctly
+    gmx_rmpbc_done(gpbc)
     return x
+
+