Commit f50aaa49 authored by Berk Onat's avatar Berk Onat
Browse files

Now pymolfile read coordinates

parent aed8364e
/* Hey emacs this is -*- C -*- and this is my editor vim.
*
* molfile.c : C and Fortran interfaces for molfile_plugins
* Copyright (c) Berk Onat <b.onat@warwick.ac.uk> 2017
*
* This program is under BSD LICENSE
*/
/*
* The code is written following the plugin test
* context of f77_molfile.c by Axel Kohlmeyer and
* in molfile_plugin/src/f77 and catdcd.c by
* Justin Gullingsrud of VMD plugins.
*/
/* Get HAVE_CONFIG_H */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <ctype.h>
/* get fixed-width types if we are using ANSI C99 */
#ifdef HAVE_STDINT_H
# include <stdint.h>
#elif (defined HAVE_INTTYPES_H)
# include <inttypes.h>
#endif
#include "pymolfile.h"
int numplugins=0;
molfile_plugin_t** plugin_list;
/* * * * * * * * * * * * * * * * * * * * * * *
* Helper functions to set and store plugins *
* * * * * * * * * * * * * * * * * * * * * * */
#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
#else
#define PyBytes_FromString PyString_FromString
#endif
#define NPY_NEXT_ALIGNED_OFFSET(offset, alignment) \
(((offset) + (alignment) - 1) & (-(alignment)))
#if 0
static NPY_INLINE npy_bool
_is_basic_python_type(PyTypeObject *tp)
{
return (
/* Basic number types */
tp == &PyBool_Type ||
#if PY_VERSION_HEX >= 0x03000000
tp == &PyInt_Type ||
#endif
tp == &PyLong_Type ||
tp == &PyFloat_Type ||
tp == &PyComplex_Type ||
/* Basic sequence types */
tp == &PyList_Type ||
tp == &PyTuple_Type ||
tp == &PyDict_Type ||
tp == &PySet_Type ||
tp == &PyFrozenSet_Type ||
tp == &PyUnicode_Type ||
tp == &PyBytes_Type ||
#if PY_VERSION_HEX >= 0x03000000
tp == &PyString_Type ||
#endif
/* other builtins */
tp == &PySlice_Type ||
tp == Py_TYPE(Py_None) ||
tp == Py_TYPE(Py_Ellipsis) ||
tp == Py_TYPE(Py_NotImplemented) ||
/* TODO: ndarray, but we can't see PyArray_Type here */
/* sentinel to swallow trailing || */
NPY_FALSE
);
}
static NPY_INLINE PyObject *
maybe_get_attr(PyObject *obj, char *name)
{
PyTypeObject *tp = Py_TYPE(obj);
PyObject *res = (PyObject *)NULL;
/* Attribute referenced by (char *)name */
if (tp->tp_getattr != NULL) {
res = (*tp->tp_getattr)(obj, name);
if (res == NULL) {
PyErr_Clear();
}
}
/* Attribute referenced by (PyObject *)name */
else if (tp->tp_getattro != NULL) {
#if PY_VERSION_HEX >= 0x03000000
PyObject *w = PyUnicode_InternFromString(name);
#else
PyObject *w = PyString_InternFromString(name);
#endif
if (w == NULL) {
return (PyObject *)NULL;
}
res = (*tp->tp_getattro)(obj, w);
Py_DECREF(w);
if (res == NULL) {
PyErr_Clear();
}
}
return res;
}
static NPY_INLINE PyObject *
maybe_set_attr(PyObject *obj, char *name, PyObject *attr)
{
PyTypeObject *tp = Py_TYPE(obj);
int status;
/* Attribute referenced by (char *)name */
if (tp->tp_setattr != NULL) {
status = (*tp->tp_setattr)(obj, name, attr);
if (status < 0) {
PyErr_Clear();
PyErr_Format(PyExc_IOError, "Error: no __array_interface__ attribute in object.");
return NULL;
}
}
/* Attribute referenced by (PyObject *)name */
else if (tp->tp_setattro != NULL) {
#if PY_VERSION_HEX >= 0x03000000
PyObject *w = PyUnicode_InternFromString(name);
#else
PyObject *w = PyString_InternFromString(name);
#endif
if (w == NULL) {
printf("Name is NULL\n");
PyErr_Clear();
PyErr_Format(PyExc_IOError, "Error: no __array_interface__ attribute in object.");
return NULL;
}
//status = (*tp->tp_setattro)(obj, w, attr);
//status = PyObject_GenericSetAttr(obj, w, attr);
status = PyObject_SetAttrString(obj, name, attr);
if (status < 0) {
PyErr_Clear();
PyErr_Format(PyExc_IOError, "Error: can not set interface attribute.");
return NULL;
}
}
return obj;
}
static int get_array_struct(PyObject* obj, int *maxndim, npy_intp *d_shape, PyArrayInterface* inter)
{
PyTypeObject *tp = Py_TYPE(obj);
PyObject* ret;
char* str;
int i;
str = "__array_struct__";
/* We do not need to check for special attributes on trivial types */
if (_is_basic_python_type(tp)) {
ret = NULL;
}
/* obj has the __array_struct__ interface */
ret = maybe_get_attr(obj, str);
if (ret != NULL) {
int nd = -1;
if (PyCapsule_CheckExact(ret)) {
PyArrayInterface *inter;
inter = (PyArrayInterface *)NpyCapsule_AsVoidPtr(ret);
if (inter->two == 2) {
nd = inter->nd;
if (nd >= 0) {
if (nd < *maxndim) {
*maxndim = nd;
}
for (i=0; i<*maxndim; i++) {
d_shape[i] = inter->shape[i];
}
}
}
}
if(nd >= 0){
return 0;
} else {
return -1;
}
} else {
return -1;
}
}
PyObject* get_array_attr_interface(PyObject* obj, char* str, int *maxndim, npy_intp *d_shape){
PyTypeObject *tp = Py_TYPE(obj);
char* dictstr;
int i, failflag=0;
dictstr = "__array_interface__";
PyObject* ret;
PyObject* retuni = NULL;
PyObject* newi = NULL;
PyObject* objj = NULL;
if (_is_basic_python_type(tp)) {
ret = NULL;
}
/* obj has the __array_interface__ interface */
ret = maybe_get_attr(obj, dictstr);
if (ret != NULL) {
int nd = -1;
if (PyDict_Check(ret)) {
objj = PyDict_GetItemString(ret, str);
if (str == "shape"){
if (objj && PyTuple_Check(objj)) {
nd = PyTuple_GET_SIZE(objj);
if (nd >= 0) {
maxndim = &nd;
for (i=0; i<nd; i++) {
retuni = PyTuple_GET_ITEM(objj, i);
}
}
}
} else {
PyObject* tuple = PyTuple_New(1);
PyTuple_SET_ITEM(tuple, 0, objj);
return (PyObject*) tuple;
}
}
if (nd >= 0) {
PyObject* tuple = PyTuple_New(1);
PyTuple_SET_ITEM(tuple, 0, objj);
return (PyObject*) tuple;
} else {
failflag = 1;
}
} else {
failflag = 1;
}
if (failflag){
PyErr_Clear();
PyErr_Format(PyExc_IOError, "Error: can not get attribute of given numpy array.");
return NULL;
}
}
PyObject* get_descr_array_interface(PyObject* obj){
PyTypeObject *tp = Py_TYPE(obj);
char* dictstr;
char* descrstr;
int i, failflag=0;
dictstr = "__array_interface__";
descrstr = "descr";
PyObject* ret;
PyObject* retuni = NULL;
PyObject* newi = NULL;
PyObject* objj = NULL;
if (_is_basic_python_type(tp)) {
ret = NULL;
}
/* obj has the __array_interface__ interface */
ret = maybe_get_attr(obj, dictstr);
if (ret != NULL) {
if (PyDict_Check(ret)) {
objj = PyDict_GetItemString(ret, descrstr);
return (PyObject*) objj;
} else {
failflag = 1;
}
} else {
failflag = 1;
}
if (failflag){
PyErr_Clear();
PyErr_Format(PyExc_IOError, "Error: can not get attribute of given numpy array.");
return NULL;
}
}
#endif
molfile_plugin_t* get_plugin(molfile_plugin_t** plug_list, int plug_no)
{
molfile_plugin_t* plugin;
if(plug_no < 0){
plugin = NULL;
} else {
plugin = plug_list[plug_no];
}
return plugin;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * *
* Interface functions to initialize molfile plugins *
* * * * * * * * * * * * * * * * * * * * * * * * * * */
/* check validity of plugins and register them. */
static int molfile_register(void*, vmdplugin_t *plugin) {
if (!plugin->type || !plugin->name || !plugin->author) {
// skipping plugin with incomplete header
return VMDPLUGIN_ERROR;
}
else if (plugin->abiversion != vmdplugin_ABIVERSION) {
// skipping plugin with incompatible ABI
return VMDPLUGIN_ERROR;
}
else if (0 != strncmp(plugin->type, "mol file", 8)) {
// skipping plugin of incompatible type
return VMDPLUGIN_ERROR;
}
else if (numplugins >= MAXPLUGINS) {
// too many plugins: increase MAXPLUGINS
return VMDPLUGIN_ERROR;
}
//if (plugin_find(&plugindict, plugin->name) != NULL) {
// multiple plugins for file type
// return VMDPLUGIN_ERROR;
//} else {
plugin_list[numplugins] = (molfile_plugin_t *) plugin;
//plugin_add(&plugindict, plugin->name, numplugins);
++numplugins;
return VMDPLUGIN_SUCCESS;
//}
}
molfile_plugin_t** molfile_plugin_list(int maxsize)
{
if(maxsize < MAXPLUGINS){
maxsize = MAXPLUGINS;
}
plugin_list = (molfile_plugin_t**) malloc(sizeof(molfile_plugin_t*)*maxsize);
return plugin_list;
}
/* register all available plugins and clear handles. */
int molfile_init(void)
{
MOLFILE_INIT_ALL;
MOLFILE_REGISTER_ALL(NULL,molfile_register);
return numplugins;
}
/* unregister all available plugins */
int molfile_finish(void)
{
MOLFILE_FINI_ALL;
return 0;
}
/* * * * * * * * * * * * * * * * * * * * * * *
* Wrappers to directly access molfile plugin*
* functions and settings *
* * * * * * * * * * * * * * * * * * * * * * */
/* molfile_plugin_t access */
/* Functions in molfile_plugin_t */
PyObject* read_fill_structure(PyObject* molpack, PyObject* prototype)
{
// Py_Initialize();
import_array();
int options = 0;
molfile_plugin_t* plugin;
void* file_handle;
molfile_atom_t* data;
int numatoms, status;
int nd;
PyObject *ret = NULL;
// Access plugin_handle values
MolObject* plugin_handle = (MolObject*) molpack;
plugin = plugin_handle->plugin;
file_handle = plugin_handle->file_handle;
numatoms = plugin_handle->natoms;
// Allocate memory for array of molfile_atom_t struct
data = (molfile_atom_t *)calloc(numatoms,sizeof(molfile_atom_t));
// Get array values in molfile_atom_t
if (plugin->read_structure) {
status = plugin->read_structure(file_handle, &options, data);
// Check if the plugin returns the results
if (status!=0){
PyErr_Format(PyExc_IOError, "Error accessing molfile_atom_t in read_structure function of plugin.");
return NULL;
}
nd = 1;
npy_intp dims[1] = { numatoms };
npy_intp strides[1] = { sizeof(molfile_atom_t) };
Py_INCREF(prototype);
ret = PyArray_NewFromDescr(Py_TYPE(prototype), PyArray_DESCR((PyArrayObject*)prototype),
nd, dims,
strides, data,
PyArray_FLAGS((PyArrayObject*)prototype), prototype);
Py_DECREF(prototype);
return (PyObject*) ret;
} else {
PyErr_Format(PyExc_AttributeError, "molfile plugin does not have read_structure function.");
Py_RETURN_NONE;
}
}
PyObject* read_fill_bonds(PyObject* molpack)
{
import_array();
int options = 0;
molfile_plugin_t* plugin;
void* file_handle;
molfile_atom_t* data;
int numatoms, status;
int nd;
PyObject *ret = NULL;
PyObject *from_arr = NULL;
PyObject *to_arr = NULL;
PyObject *bondorder_arr = NULL;
PyObject *bondtype_arr = NULL;
PyObject *bondtypename_arr = NULL;
PyArrayInterface *inter = NULL;
// Access plugin_handle values
MolObject* plugin_handle = (MolObject*) molpack;
plugin = plugin_handle->plugin;
file_handle = plugin_handle->file_handle;
numatoms = plugin_handle->natoms;
if (plugin->read_bonds) {
int nbonds, *from, *to, *bondtype, nbondtypes;
float *bondorder;
char **bondtypename;
if ((status = plugin->read_bonds(file_handle, &nbonds, &from, &to,
&bondorder, &bondtype, &nbondtypes, &bondtypename))) {
PyErr_Format(PyExc_IOError, "Error accessing read_bonds function of plugin.");
return NULL;
}
inter = (PyArrayInterface*)malloc(sizeof(PyArrayInterface));
if (inter==NULL)
return PyErr_NoMemory();
inter->flags = NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE;
ret = PyDict_New();
nd = 1;
if (nbonds>0) {
npy_intp dims[1] = { nbonds };
npy_intp strides[1] = { NPY_SIZEOF_INT };
from_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT),
nd, dims,
strides, from,
inter->flags, NULL);
PyDict_SetItemString(ret, "from", from_arr);
to_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT),
nd, dims,
strides, to,
inter->flags, NULL);
PyDict_SetItemString(ret, "to", to_arr);
if (!bondorder) {
strides[0] = { NPY_SIZEOF_FLOAT };
bondorder_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_FLOAT),
nd, dims,
strides, bondorder,
inter->flags, NULL);
PyDict_SetItemString(ret, "bondorder", bondorder_arr);
}
if (!bondtype) {
strides[0] = { NPY_SIZEOF_INT };
bondtype_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT),
nd, dims,
strides, bondtype,
inter->flags, NULL);
PyDict_SetItemString(ret, "bondtype", bondtype_arr);
}
if (!bondtypename) {
dims[0] = { nbondtypes };
strides[0] = { sizeof(NPY_STRING) };
bondtypename_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_STRING),
nd, dims,
strides, bondtypename,
inter->flags, NULL);
PyDict_SetItemString(ret, "bondtypename", bondtypename_arr);
}
return (PyObject*) ret;
} else {
Py_RETURN_NONE;
}
} else {
PyErr_Format(PyExc_AttributeError, "molfile plugin does not have read_bonds function.");
Py_RETURN_NONE;
}
}
This diff is collapsed.
......@@ -233,6 +233,9 @@ extern PyObject * read_fill_bonds(PyObject* molpack);
%feature("autodoc", "0") read_fill_angles;
extern PyObject * read_fill_angles(PyObject* molpack);
%feature("autodoc", "0") read_fill_next_timestep;
extern PyObject * read_fill_next_timestep(PyObject* molpack);
/*
#define DIM 3
typedef int ivec[DIM];
......
......@@ -169,6 +169,10 @@ def read_fill_bonds(molpack: 'PyObject *') -> "PyObject *":
def read_fill_angles(molpack: 'PyObject *') -> "PyObject *":
"""read_fill_angles(molpack) -> PyObject *"""
return _libpymolfile.read_fill_angles(molpack)
def read_fill_next_timestep(molpack: 'PyObject *') -> "PyObject *":
"""read_fill_next_timestep(molpack) -> PyObject *"""
return _libpymolfile.read_fill_next_timestep(molpack)
# This file is compatible with both classic and new-style classes.
......@@ -3834,6 +3834,22 @@ fail:
}
SWIGINTERN PyObject *_wrap_read_fill_next_timestep(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
PyObject *resultobj = 0;
PyObject *arg1 = (PyObject *) 0 ;
PyObject * obj0 = 0 ;
PyObject *result = 0 ;
if (!PyArg_ParseTuple(args,(char *)"O:read_fill_next_timestep",&obj0)) SWIG_fail;
arg1 = obj0;
result = (PyObject *)read_fill_next_timestep(arg1);
resultobj = result;
return resultobj;
fail:
return NULL;
}
static PyMethodDef SwigMethods[] = {
{ (char *)"SWIG_PyInstanceMethod_New", (PyCFunction)SWIG_PyInstanceMethod_New, METH_O, NULL},
{ (char *)"del_plugin", _wrap_del_plugin, METH_VARARGS, NULL},
......@@ -3850,6 +3866,7 @@ static PyMethodDef SwigMethods[] = {
{ (char *)"read_fill_structure", _wrap_read_fill_structure, METH_VARARGS, (char *)"read_fill_structure(molpack, prototype) -> PyObject *"},
{ (char *)"read_fill_bonds", _wrap_read_fill_bonds, METH_VARARGS, (char *)"read_fill_bonds(molpack) -> PyObject *"},
{ (char *)"read_fill_angles", _wrap_read_fill_angles, METH_VARARGS, (char *)"read_fill_angles(molpack) -> PyObject *"},
{ (char *)"read_fill_next_timestep", _wrap_read_fill_next_timestep, METH_VARARGS, (char *)"read_fill_next_timestep(molpack) -> PyObject *"},
{ NULL, NULL, 0, NULL }
};
......
......@@ -387,23 +387,23 @@ PyObject* read_fill_structure(PyObject* molpack, PyObject* prototype)
PyObject *ret = NULL;
// Access plugin_handle values
MolObject* plugin_handle = (MolObject*) molpack;
if (!plugin_handle->plugin) {
if (plugin_handle->plugin) {
plugin = plugin_handle->plugin;
} else {
PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
return NULL
return NULL;
}
if (!plugin_handle->file_handle) {
if (plugin_handle->file_handle) {
file_handle = plugin_handle->file_handle;
} else {
PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
return NULL
return NULL;
}
if (!plugin_handle->natoms) {
if (plugin_handle->natoms) {
numatoms = plugin_handle->natoms;
} else {
PyErr_Format(PyExc_IOError, "no assigned number of atoms in molfile plugin handle.");
return NULL
return NULL;
}
// Allocate memory for array of molfile_atom_t struct
data = (molfile_atom_t *)calloc(numatoms,sizeof(molfile_atom_t));
......@@ -699,3 +699,88 @@ PyObject* read_fill_angles(PyObject* molpack)
}
PyObject* read_fill_next_timestep(PyObject* molpack)
{
import_array();
molfile_timestep_t timestep;
molfile_plugin_t* plugin;
void* file_handle;
int numatoms, status;
int nd;
PyObject *ret = NULL;
// Access plugin_handle values
MolObject* plugin_handle = (MolObject*) molpack;
if (plugin_handle->plugin) {
plugin = plugin_handle->plugin;
} else {
PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
return NULL;
}
if (plugin_handle->file_handle) {
file_handle = plugin_handle->file_handle;
} else {
PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
return NULL;
}
if (plugin_handle->natoms) {
numatoms = plugin_handle->natoms;
} else {
PyErr_Format(PyExc_IOError, "no assigned number of atoms in molfile plugin handle.");
return NULL;
}
if (plugin->read_next_timestep) {
PyArrayInterface *inter = NULL;
inter = (PyArrayInterface*)malloc(sizeof(PyArrayInterface));
if (inter==NULL)
return PyErr_NoMemory();
inter->flags = NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE;
ret = PyDict_New();
timestep.coords = (float *)malloc(3*numatoms*sizeof(float));
timestep.velocities = (float *)malloc(3*numatoms*sizeof(float));
if (!(status = plugin->read_next_timestep(file_handle, numatoms, &timestep))) {
if (status == MOLFILE_EOF) {
Py_RETURN_NONE;
}
else if (status != MOLFILE_SUCCESS) {
PyErr_Format(PyExc_IOError, "Failed in calling read_next_timestep function of plugin.");
Py_RETURN_NONE;
}
else {