Skip to content

Commit

Permalink
changed xtc.Trajectory class to use utilize own xdrlib and keep old s…
Browse files Browse the repository at this point in the history
…cripts working
  • Loading branch information
dseeliger committed Aug 18, 2015
1 parent 4faf52a commit a624dc6
Show file tree
Hide file tree
Showing 5 changed files with 70 additions and 126 deletions.
1 change: 1 addition & 0 deletions .gitignore
Expand Up @@ -2,3 +2,4 @@
*~
build.sh
.localhistory
build
2 changes: 1 addition & 1 deletion pmx/__init__.py
Expand Up @@ -37,7 +37,7 @@
functionality. Take a look at the example scripts.
"""
__version__ = '1.0.0'
__version__ = '1.1.0dev'

PMX_VERSION = __version__

Expand Down
57 changes: 48 additions & 9 deletions pmx/xdrfile.py
Expand Up @@ -27,20 +27,17 @@
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
#
# Adapted by Daniel Seeliger for use in the pmx package (Aug 2015)
#

from ctypes import *
import os.path

mTrr,mNumPy=1,2
try:
from numpy import *
from numpy.ctypeslib import ndpointer
auto_mode=mNumPy
except:
auto_mode=0
auto_mode=0


class frame:
class Frame:
#variables
#x: rvec*natoms / numpy array if installed
#box DIM*DIM
Expand All @@ -51,6 +48,7 @@ class frame:

def __init__(self,n,mode):
#create vector for x
self.natoms = n
if mode&mNumPy:
self.x=empty((n,3),dtype=float32)
self.box = empty((3,3),float32)
Expand All @@ -59,7 +57,48 @@ def __init__(self,n,mode):
self.box = (c_float*3*3)()


class xdrfile:
def update_box( self, box ):
for i in range(3):
for k in range(3):
box[i][k] = self.box[i][k]

def update_atoms( self, atom_sel ):
for i, atom in enumerate(atom_sel.atoms):
if atom_sel.unity == 'A':
atom.x[0] = self.x[i][0]*10
atom.x[1] = self.x[i][1]*10
atom.x[2] = self.x[i][2]*10
else:
atom.x[0] = self.x[i][0]
atom.x[1] = self.x[i][1]
atom.x[2] = self.x[i][2]

def update( self, atom_sel ):
self.update_atoms(atom_sel )
self.update_box( atom_sel.box )


def get_natoms(self):
return self.natoms
def get_time(self):
return self.time
def get_step(self):
return self.step
def get_prec(self):
return self.prec
def get_box(self):
return self.box


def __str__(self):

s = '< xdrlib.Frame: natoms = %d | step = %d | time = %8.2f >' % (self.get_natoms(), self.get_step(), self.get_time())
return s




class XDRFile:
exdrOK, exdrHEADER, exdrSTRING, exdrDOUBLE, exdrINT, exdrFLOAT, exdrUINT, exdr3DX, exdrCLOSE, exdrMAGIC, exdrNOMEM, exdrENDOFFILE, exdrNR = range(13)

#
Expand Down Expand Up @@ -118,7 +157,7 @@ def __init__(self,fn,mode="Auto",ft="Auto"):


def __iter__(self):
f = frame(self.natoms,self.mode)
f = Frame(self.natoms,self.mode)
#temporary c_type variables (frame variables are python type)
step = c_int()
time = c_float()
Expand Down
116 changes: 6 additions & 110 deletions pmx/xtc.py
Expand Up @@ -28,125 +28,21 @@
# CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
# ----------------------------------------------------------------------

import sys, os
import sys, os, xdrfile

from ctypes import cdll, c_int, c_float, c_double
from ctypes import c_char_p,POINTER,c_int,byref

c_real = c_float
rvec=c_real*3
matrix=c_real*3*3
class Trajectory(xdrfile.XDRFile):
""" changed Trajectory class to use new XDRFile class and keep old scripts working"""

def __init__(self, filename, **kwargs):
xdrfile.XDRFile.__init__(self, filename, **kwargs)

class Trajectory:

def __init__(self, filename):
self.filename = filename
self.libgmx_path = self.__check_gmxlib()
self.libgmx = cdll.LoadLibrary(self.libgmx_path)
self.fp = None
self.natoms = c_int()
self.step = c_int()
self.time = c_real()
self.prec = c_real()
self.bOK = c_int()
self.box = matrix()
self.x = POINTER(rvec)()
self.__have_first_frame = False
self.open_xtc(filename)
self.read_first_xtc()

def __check_gmxlib(self):
p = os.environ.get('GMX_DLL')
if not p:
print >> sys.stderr, "pmx_Error_> Path to Gromacs shared libraries is not set (GMX_DLL)"
print >> sys.stderr, "pmx_Error_> Cannot load \"libgmx.so\""
sys.exit(1)
return os.path.join(p,'libgmx.so')


def open_xtc(self, filename):
self.libgmx.open_xtc.restype = POINTER(c_char_p)
self.fp = self.libgmx.open_xtc(filename, "r")

def close_xtc(self):
self.libgmx.close_xtc(self.fp)

def read_first_xtc(self):

ret = self.libgmx.read_first_xtc(self.fp,byref(self.natoms),byref(self.step),byref(self.time),self.box,byref(self.x),byref(self.prec),byref(self.bOK))
if not ret:
print >> sys.stderr, "pmx_Error_> in Trajectory.read_first_xtc()"
sys.exit(1)
self.__have_first_frame = True

def read_next_xtc(self):
if not self.__have_first_frame:
print >> sys.stderr, "pmx_Error_> First frame not read"
sys.exit(1)
ret = self.libgmx.read_next_xtc(self.fp,self.natoms,byref(self.step),byref(self.time),self.box,self.x,byref(self.prec),byref(self.bOK))
return ret

def update_box( self, box ):
for i in range(3):
for k in range(3):
box[i][k] = self.box[i][k]

def update_atoms( self, atom_sel ):
assert len(atom_sel.atoms) == self.natoms.value
for i, atom in enumerate(atom_sel.atoms):
if atom_sel.unity == 'A':
atom.x[0] = self.x[i][0]*10
atom.x[1] = self.x[i][1]*10
atom.x[2] = self.x[i][2]*10
else:
atom.x[0] = self.x[i][0]
atom.x[1] = self.x[i][1]
atom.x[2] = self.x[i][2]

def update( self, atom_sel ):
self.update_atoms(atom_sel )
self.update_box( atom_sel.box )

def get_box(self):
return [
[self.box[0][0], self.box[0][1], self.box[0][1]],
[self.box[1][0], self.box[1][1], self.box[1][1]],
[self.box[2][0], self.box[2][1], self.box[2][1]],
]

def get_natoms(self):
return self.natoms.value
def get_time(self):
return self.time.value
def get_step(self):
return self.step.value
def get_prec(self):
return self.prec.value
def get_x(self, unity='nm'):
x = []
for i in range(self.get_natoms):
if unity == 'A':
x.append( [self.x[i][0]*10, self.x[i][1]*10, self.x[i][2]*10 ])
else:
x.append( [self.x[i][0], self.x[i][1], self.x[i][2] ])
return x

def __str__(self):
return self.natoms

s = '< Trajectory: %s | natoms = %d | step = %d | time = %8.2f >' % (self.filename, self.get_natoms(), self.get_step(), self.get_time())
return s


def __iter__(self):
return self

def next(self):

ok = self.read_next_xtc()
if not ok:
raise StopIteration
return self



20 changes: 14 additions & 6 deletions setup.py
Expand Up @@ -42,9 +42,17 @@ def run(self):

pmx = Extension('pmx/_pmx',
libraries = ['m'],
include_dirs = ['src'],
sources = ['src/Geometry.c','src/wrap_Geometry.c',
'src/init.c', 'src/Energy.c']
include_dirs = ['src/pmx'],
sources = ['src/pmx/Geometry.c','src/pmx/wrap_Geometry.c',
'src/pmx/init.c', 'src/pmx/Energy.c']

)

xdrio = Extension('pmx/_xdrio',
libraries = ['m'],
include_dirs = ['src/xdr'],
sources = ['src/xdr/xdrfile.c','src/xdr/xdrfile_trr.c',
'src/xdr/xdrfile_xtc.c']

)

Expand All @@ -58,11 +66,11 @@ def run(self):


setup (name = 'pmx',
version = '1.0.0',
version = '1.1.0dev',
description = 'Python Toolbox structure file editing and writing simulation setup/analysis tools',
author = 'Daniel Seeliger',
author_email = 'seeliger.biosoft@gmail.de',
url = 'http://code.google.com/p/pmx/',
url = 'https://github.com/dseeliger/pmx/',
long_description = '''Tools to play with structure files, topologies, index files, xtc files, etc....''',
packages = ['pmx'],
data_files = [('pmx/data',['data/bbdep.pkl']),
Expand All @@ -73,7 +81,7 @@ def run(self):
('pmx/data',['data/ffamber99sbnb.itp']),
('pmx/data',['data/blosum62_new.mat'])
],
ext_modules = [pmx],
ext_modules = [pmx,xdrio],
cmdclass = {'install_data': install_data_files
},

Expand Down

0 comments on commit a624dc6

Please sign in to comment.