Commit 954ddf92 authored by Berk Onat's avatar Berk Onat

OpenMolfile Class is added with trajectory and topology functions

parent 59b6b410
from __future__ import absolute_import
from __future__ import print_function
import os
import re
import sys
import numpy as np
import warnings
if sys.version_info > (3,):
long = int
try:
from .molfile import libpymolfile
except ImportError:
warnings.warn("libpymolfile package not available, pymolfile does not work without its library!")
def plugins():
""" Information on the available molfile plugins
Example tuple: ('psf', 'psf', 1, 1, 1, 0, 1, 1, 1, 0,
'CHARMM,NAMD,XPLOR PSF', 'mol file reader',
'Justin Gullingsrud, John Stone', 1, 9, 17, 1)
The fields in the tuple represent info in ordered as follows:
1: format extension
2: format name
3: read_structure is avaliable if 1
4: read_bonds is avaliable if 1
5: read_angles is avaliable if 1
6: read_next_timestep is avaliable if 1
7: write_structure is avaliable if 1
8: write_bonds is avaliable if 1
9: write_angles is avaliable if 1
10: write_timestep is avaliable if 1
11: long name of the plugin
12: type of plugin
13: authors of the plugin
14: major version of the plugin
15: minor version of the plugin
16: ABI version of the plugin
17: 1 if is reentrant (returns is_reentrant)
Returns: A list of tuples that includes the information and
capabilities of each molfile plugin. The information is
extracted from molfile_plugin_t.
"""
max_num_plugins = 200
c_list = libpymolfile.molfile_plugin_list(max_num_plugins)
numlist = libpymolfile.molfile_init()
plugins_list = [libpymolfile.molfile_plugin_info(c_list, i) for i in range(numlist)]
libpymolfile.molfile_finish()
return plugins_list
def open(file_name_with_path, topology=None,
file_format=None, topology_format=None,
plugin=None, topology_plugin=None):
""" The main function to read topology and
trajectory files
Returns: Depends on the file format and arguments:
If structure file is supplied, it returns topology class.
If only trajectory file is supplied, it returns trajectory class
without topology information. (the number of atoms must be known)
If both files are supplied, it returns trajectory class with
topology information.
"""
pass
def get_dir_base_extension(file_name):
""" Splits directory, file base and file extensions
Returns: directory without leading '/',
file base name, and file extension without '.'
"""
file_base, file_extension_with_dot = os.path.splitext(os.path.basename(file_name))
file_extension = file_extension_with_dot.split(".")[-1]
file_dir = os.path.dirname(file_name)
return file_dir, file_base, file_extension
def get_extension(file_name):
""" Gets file extension of a file
Returns: file extension without '.'
"""
file_extension_with_dot = os.path.splitext(os.path.basename(file_name))[1]
return file_extension_with_dot.split(".")[-1]
def get_plugin_with_ext(file_ext):
""" Search molfile plugins list and returns the plugin info
for the first matched extension.
Returns: Plugin no in the list and the list item (the plugin info tuple)
"""
if not MOLFILE_PLUGINS:
MOLFILE_PLUGINS = plugins()
if MOLFILE_PLUGINS:
plugin_no = -1
for plugin_info in MOLFILE_PLUGINS:
plugin_no += 1
if file_ext == plugin_info[1]:
return (plugin_no, plugin_info)
return None
def get_plugin_with_name(plugin_name):
""" Search molfile plugins list and returns the plugin info
for the first matching name in plugin name field.
Returns: Plugin no in the list and the list item (the plugin info tuple)
"""
if not MOLFILE_PLUGINS:
MOLFILE_PLUGINS = plugins()
if MOLFILE_PLUGINS:
plugin_no = -1
for plugin_info in MOLFILE_PLUGINS:
plugin_no += 1
if plugin_name == plugin_info[2]:
return (plugin_no, plugin_info)
return None
class OpenMolfile(object):
self.plugin_list = MOLFILE_PLUGINS
self.trajectory = None
self.topology = None
self.smolobject = None
self.cmolobject = None
self.kwords = {
"file_format" : None,
"file_plugin" : None,
"topology" : None,
"topology_format" : None,
"topology_plugin" : None,
"natoms" : None
}
def __init__(self, file_name, **kwargs):
if kwargs:
for k, v in kwargs.items():
if k in self.kwords:
self.kwords[k] = v
if file_name:
if self.kwords["file_format"] is None:
file_dir, file_base, file_ext = get_dir_base_extension(file_name)
if file_ext:
self.kwords["file_format"] = file_ext
else:
self.kwords["file_format"] = file_base
if self.kwords["file_plugin"] is None:
self.kwords["file_plugin"] = "auto"
if "auto" in self.kwords["file_plugin"]:
plugin_item = get_plugin_with_ext(self.kwords["file_format"])
self.kwords["file_plugin"] =
if self.kwords["file_plugin"]:
if not MOLFILE_PLUGINS:
MOLFILE_PLUGINS = plugins()
c_list = libpymolfile.molfile_plugin_list(max_num_plugins)
numlist = libpymolfile.molfile_init()
self.cmolobject = libpymolfile.get_plugin(MOLFILE_PLUGINS, 83) #trr
else:
return None
def get_listof_parts( filename, filepath, fileformats ):
pattern = re.escape( filename[1:-4] ) + "\.part[0-9]{4,4}\.(xtc|trr)$"
parts = []
for f in os.listdir( directory ):
m = re.match( pattern, f )
if m and os.path.isfile( os.path.join( directory, f ) ):
parts.append( os.path.join( directory, f ) )
return sorted( parts )
def get_trajectory( file_name ):
ext = os.path.splitext( file_name )[1].lower()
types = {
".xtc": XtcTrajectory,
".trr": TrrTrajectory,
".netcdf": NetcdfTrajectory,
".nc": NetcdfTrajectory,
".dcd": DcdTrajectory,
}
if ext in types:
return types[ ext ]( file_name )
else:
raise Exception( "extension '%s' not supported" % ext )
class Topology( object ):
self.structure = None
self.bonds = None
self.angles = None
def __init__( self, file_name ):
pass
class Trajectory( object ):
def __init__( self, file_name ):
pass
def update( self, force=False ):
pass
def _get_frame( self, index ):
pass
def get_frame( self, index, atom_indices=None ):
box, coords, time = self._get_frame( int( index ) )
if atom_indices:
coords = np.concatenate([
])
return {
"coords": coords
"velocities": coords
}
def get_frame_string( self, index, atom_indices=None ):
frame = self.get_frame( index, atom_indices=atom_indices )
return (
)
def __del__( self ):
pass
......@@ -5,5 +5,7 @@ from . import pymolfile
__all__ = [ "pymolfile" ]
MAX_NUM_PLUGINS = 200
C_MOLFILE_PLUGINS = libpymolfile.molfile_plugin_list(MAX_NUM_PLUGINS)
MOLFILE_PLUGINS = pymolfile.plugins()
echo "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"
swig -py3 -Wall -c++ -python libpymolfile.i
echo "'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'-._.-'"
g++-7 -fPIC -Wall -Wextra -shared -g -Wunused-function -Wunused-parameter \
-I/labEnv3/lib/python3.6/site-packages/numpy/core/include/ \
-I./include \
......
......@@ -30,9 +30,6 @@
#include <ctype.h>
#include <string.h>
#include <inttypes.h>
#include "molfile_plugin.h"
#include "libmolfile_plugin.h"
#include "vmdplugin.h"
#include "pymolfile.h"
%}
......@@ -57,7 +54,7 @@ import_array();
initialize and finalize molfile plugins
*/
%feature("autodoc", "0") molfile_plugin_list;
extern molfile_plugin_t** molfile_plugin_list(int maxsize);
extern PyObject* molfile_plugin_list(int maxsize);
%feature("autodoc", "0") molfile_init;
extern int molfile_init(void);
......@@ -66,16 +63,20 @@ extern int molfile_init(void);
extern int molfile_finish(void);
%feature("autodoc", "0") get_plugin;
extern molfile_plugin_t* get_plugin(molfile_plugin_t** plugin_list, int plugin_no);
extern PyObject* get_plugin(PyObject* molcapsule, int plug_no);
%feature("autodoc", "0") molfile_plugin_info;
extern PyObject * molfile_plugin_info(PyObject* molcapsule, int plugin_no);
/*
%exception molfile_plugin_info {
$action
if (PyErr_Occurred()) SWIG_fail;
}
%inline %{
PyObject * molfile_plugin_info(molfile_plugin_t** plugin_list, int plugin_no) {
PyObject * molfile_plugin_info(PyObject* molcapsule, int plugin_no) {
molfile_plugin_t *plugin;
molfile_plugin_t** plugin_list = (molfile_plugin_t**) PyMolfileCapsule_AsVoidPtr(molcapsule);
int *plugno = &plugin_no;
int has_readstructure = 0;
int has_readbonds = 0;
......@@ -87,12 +88,12 @@ PyObject * molfile_plugin_info(molfile_plugin_t** plugin_list, int plugin_no) {
int has_writetimestep = 0;
int plugin_list_size = sizeof(plugin_list) / sizeof(molfile_plugin_t**);
if (plugno==NULL || plugin_no<0){
PyErr_Format(PyExc_IOError, "[%d] Error: molfile plugin handle no should be given, be positive value and should not exceed the list length'%d'. You set '%d'", pluginNOINIT, plugin_list_size, plugin_no);
PyErr_Format(PyExc_IOError, "Error: molfile plugin handle no should be given, be positive value and should not exceed the list length'%d'. You set '%d'", plugin_list_size, plugin_no);
return 0;
}
plugin = plugin_list[plugin_no];
if(plugin==NULL || !plugin->open_file_read){
PyErr_Format(PyExc_IOError, "[%d] Error: molfile plugin '%d' is not initialized.", pluginNOINIT, plugin_no);
PyErr_Format(PyExc_IOError, "Error: molfile plugin '%d' is not initialized.", plugin_no);
return 0;
}
if (plugin->read_structure) has_readstructure = 1;
......@@ -124,6 +125,7 @@ PyObject * molfile_plugin_info(molfile_plugin_t** plugin_list, int plugin_no) {
return tuple;
}
%}
*/
%feature("autodoc", "0") my_open_file_read;
......@@ -133,11 +135,12 @@ PyObject * molfile_plugin_info(molfile_plugin_t** plugin_list, int plugin_no) {
if (PyErr_Occurred()) SWIG_fail;
}
%inline %{
PyObject * my_open_file_read(molfile_plugin_t* plugin, char* fname, char* ftype, int natoms) {
PyObject * my_open_file_read(PyObject* molcapsule, char* fname, char* ftype, int natoms) {
if (PyType_Ready(&MolObjectType) < 0)
return NULL;
PyTypeObject *type = &MolObjectType;
MolObject *plugin_c;
molfile_plugin_t* plugin = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(molcapsule);
plugin_c = (MolObject *)type->tp_alloc(type, 0);
plugin_c->plugin = plugin;
plugin_c->file_handle = plugin->open_file_read(fname, ftype, &natoms);
......@@ -159,7 +162,7 @@ PyObject * my_close_file_read(PyObject* molpack) {
void* file_handle;
int numatoms;
MolObject* plugin_handle = (MolObject*) molpack;
plugin = plugin_handle->plugin;
plugin = plugin_handle->plugin;
file_handle = plugin_handle->file_handle;
numatoms = plugin_handle->natoms;
plugin->close_file_read(file_handle);
......
......@@ -107,8 +107,8 @@ except __builtin__.Exception:
_newclass = 0
def molfile_plugin_list(maxsize: 'int') -> "molfile_plugin_t **":
"""molfile_plugin_list(maxsize) -> molfile_plugin_t **"""
def molfile_plugin_list(maxsize: 'int') -> "PyObject *":
"""molfile_plugin_list(maxsize) -> PyObject *"""
return _libpymolfile.molfile_plugin_list(maxsize)
def molfile_init() -> "int":
......@@ -119,17 +119,17 @@ def molfile_finish() -> "int":
"""molfile_finish() -> int"""
return _libpymolfile.molfile_finish()
def get_plugin(plugin_list: 'molfile_plugin_t **', plugin_no: 'int') -> "molfile_plugin_t *":
"""get_plugin(plugin_list, plugin_no) -> molfile_plugin_t *"""
return _libpymolfile.get_plugin(plugin_list, plugin_no)
def get_plugin(molcapsule: 'PyObject *', plug_no: 'int') -> "PyObject *":
"""get_plugin(molcapsule, plug_no) -> PyObject *"""
return _libpymolfile.get_plugin(molcapsule, plug_no)
def molfile_plugin_info(plugin_list: 'molfile_plugin_t **', plugin_no: 'int') -> "PyObject *":
"""molfile_plugin_info(plugin_list, plugin_no) -> PyObject *"""
return _libpymolfile.molfile_plugin_info(plugin_list, plugin_no)
def molfile_plugin_info(molcapsule: 'PyObject *', plugin_no: 'int') -> "PyObject *":
"""molfile_plugin_info(molcapsule, plugin_no) -> PyObject *"""
return _libpymolfile.molfile_plugin_info(molcapsule, plugin_no)
def open_file_read(plugin: 'molfile_plugin_t *', fname: 'char *', ftype: 'char *', natoms: 'int') -> "PyObject *":
"""open_file_read(plugin, fname, ftype, natoms) -> PyObject *"""
return _libpymolfile.open_file_read(plugin, fname, ftype, natoms)
def open_file_read(molcapsule: 'PyObject *', fname: 'char *', ftype: 'char *', natoms: 'int') -> "PyObject *":
"""open_file_read(molcapsule, fname, ftype, natoms) -> PyObject *"""
return _libpymolfile.open_file_read(molcapsule, fname, ftype, natoms)
def close_file_read(molpack: 'PyObject *') -> "PyObject *":
"""close_file_read(molpack) -> PyObject *"""
......
This diff is collapsed.
This diff is collapsed.
......@@ -31,9 +31,6 @@ extern "C"
#include <numpy/arrayobject.h>
//enum { pluginOK, pluginNOINIT, pluginCLOSE, pluginNOMEM,
// pluginENDOFFILE, pluginFILENOTFOUND, pluginFORMATERROR };
#ifndef MAXPLUGINS
#define MAXPLUGINS 200
#endif
......@@ -46,6 +43,55 @@ struct MolObject {
MolObject(void) {}
};
static void * PyMolfileCapsule_AsVoidPtr(PyObject *obj);
static PyObject * PyMolfileCapsule_FromVoidPtr(void *ptr, void (*data)(PyObject *));
void del_molfile_plugin_list(PyObject* molcapsule);
#if PY_VERSION_HEX >= 0x03000000
#define PyInt_AsSsize_t PyLong_AsSsize_t
//#define PyInt_AsLong PyLong_AsLong
#define PyArray_Check(op) PyObject_TypeCheck(op, &PyArray_Type)
//#define PyString_FromString PyBytes_FromString
#define PyUString_Check PyUnicode_Check
#define PyUString_GET_SIZE PyUnicode_GET_SIZE
#define PyUString_FromFormat PyUnicode_FromFormat
//#define PyInt_FromLong PyLong_FromLong
#define PyString_Type PyBytes_Type
#define PyInt_Type PyLong_Type
static void * PyMolfileCapsule_AsVoidPtr(PyObject *obj)
{
void *ret = PyCapsule_GetPointer(obj, NULL);
if (ret == NULL) {
PyErr_Clear();
}
return ret;
}
static PyObject * PyMolfileCapsule_FromVoidPtr(void *ptr, void (*destr)(PyObject *))
{
PyObject *ret = PyCapsule_New(ptr, NULL, destr);
if (ret == NULL) {
PyErr_Clear();
}
return ret;
}
#else
#define PyBytes_FromString PyString_FromString
static void * PyMolfileCapsule_AsVoidPtr(PyObject *obj)
{
return PyCObject_AsVoidPtr(obj);
}
static PyObject * PyMolfileCapsule_FromVoidPtr(void *ptr, void (*destr)(void *))
{
return PyCObject_FromVoidPtr(ptr, destr);
}
#endif
static void MolObject_dealloc(MolObject* self)
{
Py_XDECREF(self->plugin);
......@@ -184,16 +230,15 @@ static PyTypeObject MolObjectType = {
MolObject_new, /* tp_new */
};
PyObject * molfile_plugin_info(PyObject* molcapsule, int plugin_no);
PyObject * read_fill_structure(PyObject* molpack, PyObject* prototype);
PyObject * read_fill_bonds(PyObject* molpack);
PyObject * read_fill_angles(PyObject* molpack);
PyObject * read_fill_next_timestep(PyObject* molpack);
PyObject * is_plugin_same(PyObject* molpack_a, PyObject* molpack_b);
PyObject * is_filehandle_same(PyObject* molpack_a, PyObject* molpack_b);
molfile_plugin_t* get_plugin(molfile_plugin_t** plug_list, int plug_no);
molfile_plugin_t** molfile_plugin_list(int maxsize);
PyObject * get_plugin(PyObject* molcapsule, int plug_no);
PyObject * molfile_plugin_list(int maxsize);
int molfile_init(void);
......
This diff is collapsed.
import numpy
import molfile.libpymolfile as pym
import pymolfile as pym
mylist = pym.molfile_plugin_list(200)
numlist = pym.molfile_init()
print(numlist)
for i in range(numlist):
testplugin = pym.molfile_plugin_info(mylist, i)
print(i, testplugin)
#splugin = pym.get_plugin(mylist, 99) #pdb
cplugin = pym.get_plugin(mylist, 83) #trr
#cplugin = pym.get_plugin(mylist, 85) #xtc
#splugin = pym.get_plugin(mylist, 105) #psf
splugin = pym.get_plugin(mylist, 81) #gro
#cplugin = pym.get_plugin(mylist, 69) #dcd
#cplugin = pym.get_plugin(mylist, 126) #nc
print(splugin)
print(cplugin)
natoms=0
#sfname="../test/DPDP.pdb"
#cfname="../test/DPDP.nc"
sfname="../test/DPDP.pdb"
cfname="../test/DPDP.nc"
#sfname="../test/ala3.pdb"
#sfname="../test/ala3.psf"
#cfname="../test/ala3.dcd"
sfname="../test/md.gro"
cfname="../test/md.trr"
#sfname="../test/md.gro"
#cfname="../test/md.trr"
#cfname="../test/md.xtc"
#sfname="../test/md_1u19.gro"
#cfname="../test/md_1u19.xtc"
#sftype="pdb"
cftype="trr"
#cftype="xtc"
#sftype="psf"
sftype="gro"
#cftype="dcd"
#cftype="nc"
print("Structure plugin:")
spluginhandle = pym.open_file_read(splugin, sfname, sftype, natoms)
print("Coordinate plugin:")
cpluginhandle = pym.open_file_read(cplugin, cfname, cftype, natoms)
print(cpluginhandle)
print(cpluginhandle.natoms)
print(spluginhandle)
print(spluginhandle.natoms)
if cpluginhandle.natoms != spluginhandle.natoms:
print("Number of atoms do not match in structure and coordinate files.")
if pym.is_plugin_same(spluginhandle, cpluginhandle):
pym.close_file_read(cpluginhandle)
cpluginhandle = spluginhandle
#if 0 && vmdplugin_ABIVERSION > 17
# /* The new PDB file formats allows for much larger structures, */
# /* which can therefore require longer chain ID strings. The */
# /* new PDBx/mmCIF file formats do not have length limits on */
# /* fields, so PDB chains could be arbitrarily long strings */
# /* in such files. At present, we know we need at least 3-char */
# /* chains for existing PDBx/mmCIF files. */
# char chain[4]; /**< required chain name, or "" */
#else
# char chain[2]; /**< required chain name, or "" */
#endif
# Change 'chain', 'S4' to S2 for PDB
if 'pdb' in [sftype, cftype]:
print("file chain is set to S2")
chain_size = 'S2'
else:
chain_size = 'S4'
x = numpy.array([
('C1','C','ACE',0,'','','','',1.0,1.0,1.0,1.0,1.0,6),
('C2','C','ACE',0,'','','','',1.0,1.0,1.0,1.0,1.0,6)
],
dtype=[
('name', 'S16'), ('type', 'S16'), ('resname', 'S8'),
('resid', 'i4'), ('segid', 'S8'), ('chain', chain_size),
('altloc', 'S2'), ('insertion', 'S2'), ('occupancy', 'f4'),
('bfactor', 'f4'), ('mass', 'f4'), ('charge', 'f4'),
('radius', 'f4'), ('atomicnumber', 'i4')
]
)
y = pym.read_fill_structure(spluginhandle, x)
print(y)
if y is not None:
print(len(y))
try:
z = pym.read_fill_bonds(spluginhandle)
print(z)
if z is not None:
print(len(z))
except (AttributeError, SystemError):
pass
try:
a = pym.read_fill_angles(spluginhandle)
print(a)
if a is not None:
print(len(a))
except (AttributeError, SystemError):
pass
if pym.is_plugin_same(spluginhandle, cpluginhandle):
pass
else:
pym.close_file_read(spluginhandle)
#step=0
#while True:
# try:
# c = pym.read_fill_next_timestep(cpluginhandle)
# if c is None:
# break
# step=step+1
# print("Step:",step)
# print(c)
# if c["coords"] is not None:
# print(len(c["coords"]))
# except (AttributeError,OSError):
# pass
pym.close_file_read(cpluginhandle)
pym.molfile_finish()
molfile = pym.OpenMolfile(sfname)
print(molfile)
import numpy
import molfile.libpymolfile as pym
mylist = pym.molfile_plugin_list(200)
print(mylist)
numlist = pym.molfile_init()
print(numlist)
for i in range(numlist):
testplugin = pym.molfile_plugin_info(mylist, i)
print(i, testplugin)
splugin = pym.get_plugin(mylist, 99) #pdb
#cplugin = pym.get_plugin(mylist, 83) #trr
#cplugin = pym.get_plugin(mylist, 85) #xtc
#splugin = pym.get_plugin(mylist, 105) #psf
#splugin = pym.get_plugin(mylist, 81) #gro
#cplugin = pym.get_plugin(mylist, 69) #dcd
cplugin = pym.get_plugin(mylist, 126) #nc
print(splugin)
print(cplugin)
natoms=0
sfname="../test/DPDP.pdb"
cfname="../test/DPDP.nc"
#sfname="../test/ala3.pdb"
#sfname="../test/ala3.psf"
#cfname="../test/ala3.dcd"
#sfname="../test/md.gro"
#cfname="../test/md.trr"
#cfname="../test/md.xtc"<