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

Added bonds, angles, dihedrals

parent f2259af4
/* 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.
......@@ -21,7 +21,7 @@
%{
/* Python SWIG interface to libpymolfile
Copyright (c) 2017 Berk Onat <b.onat@warwick.ac.uk>
Published under BSD LICENSE
Published with UIUC LICENSE
*/
#define SWIG_FILE_WITH_INIT
#define __STDC_FORMAT_MACROS
......@@ -175,6 +175,28 @@ PyObject * my_open_file_read(molfile_plugin_t* plugin, char* fname, char* ftype,
}
%}
%feature("autodoc", "0") my_close_file_read;
%rename (close_file_read) my_close_file_read;
%exception my_close_file_read {
$action
if (PyErr_Occurred()) SWIG_fail;
}
%inline %{
PyObject * my_close_file_read(PyObject* molpack) {
molfile_plugin_t* plugin;
void* file_handle;
int numatoms;
MolObject* plugin_handle = (MolObject*) molpack;
plugin = plugin_handle->plugin;
file_handle = plugin_handle->file_handle;
numatoms = plugin_handle->natoms;
plugin->close_file_read(file_handle);
Py_DECREF(plugin_handle);
Py_DECREF(molpack);
Py_RETURN_NONE;
}
%}
/*
%typemap( argout ) ( char **MolfileAtomT_CharArray )
{
......@@ -202,16 +224,22 @@ PyObject * my_open_file_read(molfile_plugin_t* plugin, char* fname, char* ftype,
}
*/
%feature("autodoc", "0") read_fill_structure;
extern PyObject * read_fill_structure(PyObject* molpack, PyObject* prototype);
%feature("autodoc", "0") read_fill_bonds;
extern PyObject * read_fill_bonds(PyObject* molpack);
%feature("autodoc", "0") read_fill_angles;
extern PyObject * read_fill_angles(PyObject* molpack);
/*
#define DIM 3
typedef int imatrix[DIM][DIM];
typedef int ivec[DIM];
typedef float fmatrix[DIM][DIM];
typedef float fvec[DIM];
*/
/* Reading from xdr files */
/*
%apply (float INPLACE_ARRAY2[ANY][ANY]) {(matrix box)}
%apply (int DIM1, int DIM2, float* INPLACE_ARRAY2) {(int natoms, int _DIM, float *x),
(int vnatoms, int v_DIM, float *v),
......@@ -224,70 +252,23 @@ typedef float fvec[DIM];
%apply (int DIM1, float* INPLACE_FARRAY2) {(int mnatoms, float* MolAtom_mass),
(int anatoms, float* MolAtom_atomicnumber)}
*/
/*
%inline %{
int structure_read(molfile_plugin_t* plugin, void* fhandle, int *options,
int* natoms, char** MolAtom_name,
int* natoms, char** MolAtom_type,
int* natoms, char** MolAtom_resname,
int* natoms, int** MolAtom_resid,
int* natoms, char** MolAtom_segid,
int* natoms, char** MolAtom_chain,
int* natoms, char** MolAtom_altloc,
int* natoms, char** MolAtom_insertion,
int* natoms, float** MolAtom_occupancy,
int* natoms, float** MolAtom_bfactor,
int* natoms, float** MolAtom_mass,
int* natoms, float** MolAtom_charge,
int* natoms, float** MolAtom_radius,
int* natoms, float** MolAtom_atomicnumber,
int* natoms, float** MolAtom_ctnumber) {
molfile_atom_t* atoms;
atoms = (molfile_atom_t *)calloc(natoms,sizeof(molfile_atom_t));
plugin->read_structure(fhandle, options, atoms);
MolAtom_name = (char **)malloc(natoms,sizeof(char*));
MolAtom_type = (char **)malloc(natoms,sizeof(char*));
MolAtom_resname = (char **)malloc(natoms,sizeof(char*));
MolAtom_resid = (int **)malloc(natoms,sizeof(int*));
%}
*/
/*
/*
%feature("autodoc", "read_bonds(MOLOBJECT, nbonds, from, to, bondorder, bondtype, nbondtypes, bondtypename) -> status") my_read_bonds;
%rename (read_bonds) my_read_bonds;
%inline %{
int structure_read(molfile_plugin_t* plugin, void* fhandle, int *options,
int natoms, char** MolAtom_name,
int tnatoms, char** MolAtom_type,
int rnatoms, int* MolAtom_resid,
int mnatoms, float* MolAtom_mass,
int anatoms, float* MolAtom_atomicnumber) {
int i;
molfile_atom_t* atoms;
molfile_atom_t atm;
atoms = (molfile_atom_t *)calloc(natoms,sizeof(molfile_atom_t));
plugin->read_structure(fhandle, options, atoms);
if(atoms == NULL) { free(atoms); return 1; }
if(atoms->type == NULL || atoms->name == NULL){ free(atoms); return 1; }
MolAtom_name = (char **)malloc(natoms,sizeof(char*));
MolAtom_type = (char **)malloc(natoms,sizeof(char*));
MolAtom_resid = (int *)malloc(natoms,sizeof(int));
MolAtom_mass = (float *)malloc(natoms,sizeof(float));
MolAtom_atomicnumber = (float *)malloc(natoms,sizeof(float));
for (i=0;i<natoms;i++){
atm = atoms[i];
MolAtom_name[i] = atm.name;
MolAtom_type[i] = atm.type;
MolAtom_resid[i] = atm.resid;
MolAtom_mass[i] = atm.mass;
MolAtom_atomicnumber[i] = atm.atomicnumber;
}
return 0;
}
PyObject * read_bonds(PyObject* molpack, int f_DIM, int *from,
int t_DIM, int *to,
int bo_DIM, float *bondorder,
int bt_DIM, int *bondtype,
int n_DIM, char *bondtypename){
(*read_bonds)(void *, int *nbonds, int **from, int **to, float **bondorder,
int **bondtype, int *nbondtypes, char ***bondtypename);
}
%}
*/