pymolfile.c 28.8 KB
Newer Older
Berk Onat's avatar
Berk Onat committed
1
2
/* Hey emacs this is -*- C -*- and this is my editor vim.
 * 
3
 * pymolfile.c : C and Python interfaces for molfile_plugins
Berk Onat's avatar
Berk Onat committed
4
5
 * Copyright (c) Berk Onat <b.onat@warwick.ac.uk> 2017
 *
6
 * This program is under UIUC Open Source License please see LICENSE file
Berk Onat's avatar
Berk Onat committed
7
8
9
10
 */

/*
 * The code is written following the plugin test 
11
 * context of main.c in molfile_plugin/src/ 
Berk Onat's avatar
Berk Onat committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
 */

/* 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"

34
35
int numplugins=0;
molfile_plugin_t** plugin_list;
36

37
38
39
40
41
#if PY_VERSION_HEX >= 0x03000000
#define PyInt_AsLong PyLong_AsLong
#define PyString_FromString PyBytes_FromString
#define PyInt_FromLong PyLong_FromLong

42
43
44
45
46
47
void del_molfile_plugin_list(PyObject* molcapsule)
{
    plugin_list = (molfile_plugin_t**) PyMolfileCapsule_AsVoidPtr(molcapsule);   
    free(plugin_list); 
    Py_XDECREF(molcapsule);
}
48
49
50
51
52
53
void del_molfile_file_handle(PyObject* molcapsule)
{
    void *file_handle = (void*) PyMolfileCapsule_AsVoidPtr(molcapsule);   
    free(file_handle); 
    Py_XDECREF(molcapsule);
}
54
#else
55
56
57
58
59
60
void* del_molfile_plugin_list(void* molcapsule)
{
    plugin_list = PyMolfileCapsule_AsVoidPtr((PyObject*)molcapsule);   
    free(plugin_list); 
    Py_XDECREF(molcapsule);
}
61
62
63
64
65
66
void* del_molfile_file_handle(void* molcapsule)
{
    void *file_handle = PyMolfileCapsule_AsVoidPtr((PyObject*)molcapsule);   
    free(file_handle); 
    Py_XDECREF(molcapsule);
}
67
68
#endif

69
70
71
/* * * * * * * * * * * * * * * * * * * * * * *
 * Helper functions to set and store plugins *
 * * * * * * * * * * * * * * * * * * * * * * */
72

73
PyObject* get_plugin(PyObject* molcapsule, int plug_no)
Berk Onat's avatar
Berk Onat committed
74
{
75
    molfile_plugin_t* plugin;
76
    molfile_plugin_t** plug_list = (molfile_plugin_t**) PyMolfileCapsule_AsVoidPtr(molcapsule);   
Berk Onat's avatar
Berk Onat committed
77
    if(plug_no < 0){
78
	plugin = NULL;
Berk Onat's avatar
Berk Onat committed
79
    } else {
80
81
82
83
84
	if(plug_list != NULL){
	    plugin = plug_list[plug_no];
	} else {
	    plugin = NULL;
	}
Berk Onat's avatar
Berk Onat committed
85
    }
86
    return PyMolfileCapsule_FromVoidPtr((void *)plugin, NULL);
Berk Onat's avatar
Berk Onat committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * *
 * 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;
    }

        plugin_list[numplugins] = (molfile_plugin_t *) plugin;
        ++numplugins;
        return VMDPLUGIN_SUCCESS;
}

117
PyObject* molfile_plugin_list(int maxsize)
Berk Onat's avatar
Berk Onat committed
118
119
120
121
122
{
    if(maxsize < MAXPLUGINS){
        maxsize = MAXPLUGINS;
    }
    plugin_list = (molfile_plugin_t**) malloc(sizeof(molfile_plugin_t*)*maxsize);
123
    return PyMolfileCapsule_FromVoidPtr((void *)plugin_list, del_molfile_plugin_list);
Berk Onat's avatar
Berk Onat committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
}

/* 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 */

150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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;
    int has_readangles = 0;
    int has_writestructure = 0;
    int has_writebonds = 0;
    int has_writeangles = 0;
    int has_readnexttimestep = 0;
    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, "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, "Error: molfile plugin '%d' is not initialized.", plugin_no);
      return 0;
    }
    if (plugin->read_structure) has_readstructure = 1;
    if (plugin->read_bonds) has_readbonds = 1;
    if (plugin->read_angles) has_readangles = 1;
    if (plugin->read_next_timestep) has_readnexttimestep = 1;
    if (plugin->write_structure) has_writestructure = 1;
    if (plugin->write_bonds) has_writebonds = 1;
    if (plugin->write_angles) has_writeangles = 1;
    if (plugin->write_timestep) has_writetimestep = 1;
    PyObject *tuple = PyTuple_New(17);
    PyTuple_SET_ITEM(tuple, 0, PyString_FromString(plugin->filename_extension));
    PyTuple_SET_ITEM(tuple, 1, PyString_FromString(plugin->name));
    PyTuple_SET_ITEM(tuple, 2, PyInt_FromLong((long)has_readstructure));
    PyTuple_SET_ITEM(tuple, 3, PyInt_FromLong((long)has_readbonds));
    PyTuple_SET_ITEM(tuple, 4, PyInt_FromLong((long)has_readangles));
    PyTuple_SET_ITEM(tuple, 5, PyInt_FromLong((long)has_readnexttimestep));
    PyTuple_SET_ITEM(tuple, 6, PyInt_FromLong((long)has_writestructure));
    PyTuple_SET_ITEM(tuple, 7, PyInt_FromLong((long)has_writebonds));
    PyTuple_SET_ITEM(tuple, 8, PyInt_FromLong((long)has_writeangles));
    PyTuple_SET_ITEM(tuple, 9, PyInt_FromLong((long)has_writetimestep));
    PyTuple_SET_ITEM(tuple, 10, PyString_FromString(plugin->prettyname));
    PyTuple_SET_ITEM(tuple, 11, PyString_FromString(plugin->type));
    PyTuple_SET_ITEM(tuple, 12, PyString_FromString(plugin->author));
    PyTuple_SET_ITEM(tuple, 13, PyInt_FromLong((long)plugin->majorv));
    PyTuple_SET_ITEM(tuple, 14, PyInt_FromLong((long)plugin->minorv));
    PyTuple_SET_ITEM(tuple, 15, PyInt_FromLong((long)plugin->abiversion));
    PyTuple_SET_ITEM(tuple, 16, PyInt_FromLong((long)plugin->is_reentrant));
    return tuple;
}

201
202
PyObject* read_fill_structure(PyObject* molpack, PyObject* prototype)
{
Berk Onat's avatar
Berk Onat committed
203
//    Py_Initialize();
204
205
206
207
208
209
210
211
212
213
    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;
Berk Onat's avatar
Berk Onat committed
214
    if (plugin_handle->plugin) {
215
216
        plugin = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle->plugin);
        //plugin = plugin_handle->plugin;   
217
    } else {
Berk Onat's avatar
Berk Onat committed
218
        PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
Berk Onat's avatar
Berk Onat committed
219
	return NULL;
Berk Onat's avatar
Berk Onat committed
220
    } 
Berk Onat's avatar
Berk Onat committed
221
    if (plugin_handle->file_handle) {
222
223
        file_handle = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle->file_handle);
        //file_handle = plugin_handle->file_handle;
224
    } else {
Berk Onat's avatar
Berk Onat committed
225
        PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
Berk Onat's avatar
Berk Onat committed
226
	return NULL; 
Berk Onat's avatar
Berk Onat committed
227
    } 
Berk Onat's avatar
Berk Onat committed
228
    if (plugin_handle->natoms) {
Berk Onat's avatar
Berk Onat committed
229
230
231
        numatoms = plugin_handle->natoms;
    } else { 
        PyErr_Format(PyExc_IOError, "no assigned number of atoms in molfile plugin handle.");
Berk Onat's avatar
Berk Onat committed
232
	return NULL;
Berk Onat's avatar
Berk Onat committed
233
    } 
234
235
236
    // 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
237
238
239
240
241
242
243
    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;
        }
Berk Onat's avatar
Berk Onat committed
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
	if(numatoms>0){
            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, "plugin read_structure does not have atoms information.");
	    Py_RETURN_NONE;
	}
259
260
261
    } else {
        PyErr_Format(PyExc_AttributeError, "molfile plugin does not have read_structure function.");
	Py_RETURN_NONE;
Berk Onat's avatar
Berk Onat committed
262
    }
263
264
}

Berk Onat's avatar
Berk Onat committed
265
266
267
268
269
270
271
272
273
274
275
276
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;
    // Access plugin_handle values
    MolObject* plugin_handle = (MolObject*) molpack;
277
    if (plugin_handle->plugin) {
278
279
        plugin = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle->plugin);
        //plugin = plugin_handle->plugin;   
280
281
282
283
284
    } else {
        PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
	return NULL;
    } 
    if (plugin_handle->file_handle) {
285
286
        file_handle = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle->file_handle);
        //file_handle = plugin_handle->file_handle;
287
288
289
290
    } else {
        PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
	return NULL; 
    } 
Berk Onat's avatar
Berk Onat committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    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;
        }
        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();
        nd = 1;
        npy_intp istrides[1] = { NPY_SIZEOF_INT };
        npy_intp fstrides[1] = { NPY_SIZEOF_FLOAT };
        npy_intp cstrides[1] = { sizeof(NPY_STRING) };
	if (nbonds>0) {
            PyObject *from_arr = NULL;
            PyObject *to_arr = NULL;
            npy_intp dims[1] = { nbonds };
            from_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                            nd, dims,
                                            istrides, from, 
                                            inter->flags, NULL);
	    PyDict_SetItemString(ret, "from", from_arr);
            to_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                          nd, dims,
                                          istrides, to, 
                                          inter->flags, NULL);
	    PyDict_SetItemString(ret, "to", to_arr);
	    if (bondorder!=NULL) {
                PyObject *bondorder_arr = NULL;
                bondorder_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_FLOAT), 
                                                     nd, dims,
                                                     fstrides, bondorder, 
			                             inter->flags, NULL);
	        PyDict_SetItemString(ret, "bondorder", bondorder_arr);
	    }
	    if (bondtype!=NULL && nbondtypes>0) {
                PyObject *bondtype_arr = NULL;
                bondtype_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                    nd, dims,
                                                    istrides, bondtype, 
			                            inter->flags, NULL);
	        PyDict_SetItemString(ret, "bondtype", bondtype_arr);
	    }
	    if (bondtypename!=NULL && nbondtypes>0) {
                PyObject *bondtypename_arr = NULL;
                npy_intp cdims[1] = { nbondtypes };
                bondtypename_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_STRING), 
                                                        nd, cdims,
                                                        cstrides, 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;
    }
}

PyObject* read_fill_angles(PyObject* molpack)
{
    import_array();
    int options = 0;
    molfile_plugin_t* plugin;
    void* file_handle;
    molfile_atom_t* data;
    int numatoms, status;
    int nd;
    int nodata = 0;
    PyObject *ret = NULL;
    // Access plugin_handle values
    MolObject* plugin_handle = (MolObject*) molpack;
373
    if (plugin_handle->plugin) {
374
375
        plugin = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle->plugin);
        //plugin = plugin_handle->plugin;   
376
377
378
379
380
    } else {
        PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
	return NULL;
    } 
    if (plugin_handle->file_handle) {
381
382
        file_handle = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle->file_handle);
        //file_handle = plugin_handle->file_handle;
383
384
385
386
    } else {
        PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
	return NULL; 
    } 
Berk Onat's avatar
Berk Onat committed
387
388
389
390
    numatoms = plugin_handle->natoms;
    // Check if there is read_angles support in this plugin
    if (plugin->read_angles) {
	// Angles
391
392
393
        int numangles;
        int *angles = NULL;
	int *angletypes = NULL;
Berk Onat's avatar
Berk Onat committed
394
        int numangletypes;
395
	char **angletypenames = NULL; 
Berk Onat's avatar
Berk Onat committed
396
	// Dihedrals
397
398
399
	int numdihedrals; 
	int *dihedrals = NULL;
	int *dihedraltypes = NULL;
Berk Onat's avatar
Berk Onat committed
400
	int numdihedraltypes;
401
        char **dihedraltypenames = NULL; 
Berk Onat's avatar
Berk Onat committed
402
	// Impropers
403
404
405
	int numimpropers;
        int *impropers = NULL;
        int *impropertypes = NULL;
Berk Onat's avatar
Berk Onat committed
406
	int numimpropertypes;
407
	char **impropertypenames = NULL;
Berk Onat's avatar
Berk Onat committed
408
	// Cterms
409
410
        int numcterms, ctermcols, ctermrows;
	int *cterms = NULL; 
Berk Onat's avatar
Berk Onat committed
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
	// Initilize zeros to number of angles, dihedrals, so on ...
	numangles = 0;
	numangletypes = 0;
	numdihedrals = 0;
	numdihedraltypes = 0;
	numimpropers = 0;
	numimpropertypes = 0;
	numcterms = 0;
	// Calling read_angles to gather the information
        if ((status = plugin->read_angles(file_handle, &numangles, &angles, &angletypes,
                                          &numangletypes, &angletypenames, &numdihedrals,
                                          &dihedrals, &dihedraltypes, &numdihedraltypes,
                                          &dihedraltypenames, &numimpropers, &impropers,        
                                          &impropertypes, &numimpropertypes, &impropertypenames,
                                          &numcterms, &cterms, &ctermcols, &ctermrows))) {
            PyErr_Format(PyExc_IOError, "Error accessing read_angles function of plugin.");
            return NULL;
        }
        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();
        nd = 1;
        npy_intp istrides[1] = { NPY_SIZEOF_INT };
        npy_intp sstrides[1] = { sizeof(NPY_STRING) };
	if (numangles>0 && angles!=NULL) {
	    nodata = 1;
            PyObject *angles_arr = NULL;
            npy_intp adims[1] = { numangles };
            angles_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                              nd, adims,
                                              istrides, angles, 
                                              inter->flags, NULL);
	    PyDict_SetItemString(ret, "angles", angles_arr);
	    if (angletypes!=NULL) {
                PyObject *angletypes_arr = NULL;
                angletypes_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                      nd, adims,
                                                      istrides, angletypes, 
                                                      inter->flags, NULL);
	        PyDict_SetItemString(ret, "angletypes", angletypes_arr);
	    }
	}
	if (numangletypes>0 && angletypenames!=NULL) {
	    nodata = 1;
            PyObject *angletypenames_arr = NULL;
            npy_intp atdims[1] = { numangletypes };
            angletypenames_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_STRING), 
                                                      nd, atdims,
                                                      sstrides, angletypenames, 
			                              inter->flags, NULL);
	    PyDict_SetItemString(ret, "angletypenames", angletypenames_arr);
	}
	if (numdihedrals>0 && dihedrals!=NULL) {
	    nodata = 1;
            PyObject *dihedrals_arr = NULL;
            npy_intp ddims[1] = { numdihedrals };
            dihedrals_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                 nd, ddims,
                                                 istrides, dihedrals, 
			                         inter->flags, NULL);
	    PyDict_SetItemString(ret, "dihedrals", dihedrals_arr);
	    if (dihedraltypes!=NULL) {
                PyObject *dihedraltypes_arr = NULL;
                dihedraltypes_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                    nd, ddims,
                                                    istrides, dihedraltypes, 
                                                    inter->flags, NULL);
	        PyDict_SetItemString(ret, "dihedraltypes", dihedraltypes_arr);
	    }
	}
	if (numdihedraltypes>0 && dihedraltypenames!=NULL) {
            PyObject *dihedraltypenames_arr = NULL;
            npy_intp dtdims[1] = { numdihedraltypes };
            dihedraltypenames_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_STRING), 
                                                         nd, dtdims,
                                                         sstrides, dihedraltypenames, 
	                                                 inter->flags, NULL);
	    PyDict_SetItemString(ret, "dihedraltypenames", dihedraltypenames_arr);
	}
	if (numimpropers>0 && impropers!=NULL) {
	    nodata = 1;
            PyObject *impropers_arr = NULL;
            npy_intp idims[1] = { numimpropers };
            impropers_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                 nd, idims,
                                                 istrides, impropers, 
			                         inter->flags, NULL);
	    PyDict_SetItemString(ret, "impropers", impropers_arr);
	    if (impropertypes!=NULL) {
                PyObject *impropertypes_arr = NULL;
                impropertypes_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                                         nd, idims,
                                                         istrides, impropertypes, 
                                                         inter->flags, NULL);
	        PyDict_SetItemString(ret, "impropertypes", impropertypes_arr);
	    }
	}
	if (numimpropertypes>0 && impropertypenames!=NULL) {
            PyObject *impropertypenames_arr = NULL;
            npy_intp itdims[1] = { numimpropertypes };
            impropertypenames_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_STRING), 
                                                         nd, itdims,
                                                         sstrides, impropertypenames, 
	                                                 inter->flags, NULL);
	    PyDict_SetItemString(ret, "impropertypenames", impropertypenames_arr);
	}
	if (numcterms>0 && cterms!=NULL) {
	    nodata = 1;
	    int ctermnd;
            npy_intp *ctermdims;
	    npy_intp *ctermstrides;
            PyObject *cterms_arr = NULL;
	    if (ctermrows>0 || ctermcols>0) {
		ctermnd = 2;
                ctermdims = (npy_intp*)calloc(ctermnd,sizeof(int));
                ctermstrides = (npy_intp*)calloc(ctermnd,sizeof(int));
	        ctermdims[0] = ctermrows;
	        ctermdims[1] = ctermcols;
	        ctermstrides[0] = NPY_SIZEOF_INT;
	        ctermstrides[1] = ctermcols*NPY_SIZEOF_INT;
	    } else {
		ctermnd = 1;
                ctermdims = (npy_intp*)calloc(ctermnd,sizeof(int));
                ctermstrides = (npy_intp*)calloc(ctermnd,sizeof(int));
	        ctermdims[0] = 8*numcterms;
	        ctermstrides[0] = NPY_SIZEOF_INT;
	    }
            cterms_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_INT), 
                                              ctermnd, ctermdims,
                                              ctermstrides, cterms, 
			                      inter->flags, NULL);
	    PyDict_SetItemString(ret, "cterms", cterms_arr);
	}
	if (nodata>0) {
            return (PyObject*) ret;
	} else {
            Py_RETURN_NONE;
	}
    } else {
        PyErr_Format(PyExc_AttributeError, "molfile plugin does not have read_angles function.");
        Py_RETURN_NONE;
    }
}

Berk Onat's avatar
Berk Onat committed
558

Berk Onat's avatar
Berk Onat committed
559
560
561
562
563
564
565
PyObject* read_fill_next_timestep(PyObject* molpack)
{
    import_array();
    molfile_plugin_t* plugin;
    void* file_handle;
    int numatoms, status;
    int nd;
566
    int i, d;
Berk Onat's avatar
Berk Onat committed
567
568
569
570
    PyObject *ret = NULL;
    // Access plugin_handle values
    MolObject* plugin_handle = (MolObject*) molpack;
    if (plugin_handle->plugin) {
571
572
        plugin = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle->plugin);
        //plugin = plugin_handle->plugin;   
573
    } else {
Berk Onat's avatar
Berk Onat committed
574
        PyErr_Format(PyExc_IOError, "molfile plugin is not active.");
575
	return NULL;
Berk Onat's avatar
Berk Onat committed
576
577
    } 
    if (plugin_handle->file_handle) {
578
579
        file_handle = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle->file_handle);
        //file_handle = plugin_handle->file_handle;
580
    } else {
Berk Onat's avatar
Berk Onat committed
581
        PyErr_Format(PyExc_IOError, "no file handle in molfile plugin handle.");
582
	return NULL; 
Berk Onat's avatar
Berk Onat committed
583
584
585
586
587
588
589
590
591
592
593
594
595
596
    } 
    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();
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
        molfile_timestep_t timestep;
	// Check if the velocities will be supplied from this plugin
	int has_velocities = -1;
	unsigned int total_steps = 0;
	unsigned int bytes_per_step = 0;
	molfile_timestep_metadata_t timestep_metadata;
        if (plugin->read_timestep_metadata) {
	    plugin->read_timestep_metadata(file_handle, &timestep_metadata);
            total_steps = timestep_metadata.count; 
	    has_velocities = timestep_metadata.has_velocities;
	    bytes_per_step = timestep_metadata.avg_bytes_per_timestep;
	} else {
	    total_steps = 0;
	    has_velocities = -2;
	}
	timestep.A=-1;
	timestep.B=-1;
	timestep.C=-1;
	timestep.alpha=-1;
	timestep.beta=-1;
	timestep.gamma=-1;
	timestep.physical_time=-1;
Berk Onat's avatar
Berk Onat committed
619
        timestep.coords = (float *)malloc(3*numatoms*sizeof(float));
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
	if(has_velocities == -2 || has_velocities == 1) { 
            timestep.velocities = (float *)malloc(3*numatoms*sizeof(float));
	    for(i=0;i<numatoms;i++)
	        for(d=0;d<3;d++)
	            timestep.velocities[i*3+d] = -1111*(d+1);
	}
        status = plugin->read_next_timestep(file_handle, numatoms, &timestep);
        if (status == MOLFILE_EOF) {
	    Py_RETURN_NONE;
	}
	else if (status != MOLFILE_SUCCESS) {
            PyErr_Format(PyExc_AttributeError, "Failed in calling read_next_timestep function of plugin.");
	    Py_RETURN_NONE;
        } 
	else {
            nd = 2;
            PyObject *coords_arr = NULL;
            npy_intp dims[2] = { numatoms, 3 };
            npy_intp strides[2] = { 3*sizeof(float), sizeof(float) };
            coords_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_FLOAT), 
Berk Onat's avatar
Berk Onat committed
640
641
642
                                              nd, dims,
                                              strides, timestep.coords, 
	                                      inter->flags, NULL);
643
644
645
646
647
	    PyDict_SetItemString(ret, "coords", coords_arr);
	    if (timestep.velocities!=NULL && (has_velocities == -2 || has_velocities == 1)) {
		if (-1111 != (int)timestep.velocities[0] && 
		    -2222 != (int)timestep.velocities[1] && 
		    -3333 != (int)timestep.velocities[2]) {
Berk Onat's avatar
Berk Onat committed
648
649
                    PyObject *velocities_arr = NULL;
                    velocities_arr = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_FLOAT), 
650
651
652
                                                          nd, dims,
                                                          strides, timestep.velocities, 
	                                                  inter->flags, NULL);
Berk Onat's avatar
Berk Onat committed
653
	            PyDict_SetItemString(ret, "velocities", velocities_arr);
654
655
656
		}
	    } else {
		if(has_velocities == -2 || has_velocities == 1)
Berk Onat's avatar
Berk Onat committed
657
		    free(timestep.velocities);
658
659
	    }
            if(timestep.A > -1)
Berk Onat's avatar
Berk Onat committed
660
	        PyDict_SetItemString(ret, "A", PyFloat_FromDouble((double)timestep.A));
661
            if(timestep.B > -1)
Berk Onat's avatar
Berk Onat committed
662
	        PyDict_SetItemString(ret, "B", PyFloat_FromDouble((double)timestep.B));
663
            if(timestep.C > -1)
Berk Onat's avatar
Berk Onat committed
664
	        PyDict_SetItemString(ret, "C", PyFloat_FromDouble((double)timestep.C));
665
            if(timestep.alpha > -1)
Berk Onat's avatar
Berk Onat committed
666
	        PyDict_SetItemString(ret, "alpha", PyFloat_FromDouble((double)timestep.alpha));
667
	    if(timestep.beta > -1)
Berk Onat's avatar
Berk Onat committed
668
	        PyDict_SetItemString(ret, "beta", PyFloat_FromDouble((double)timestep.beta));
669
	    if(timestep.gamma > -1)
Berk Onat's avatar
Berk Onat committed
670
	        PyDict_SetItemString(ret, "gamma", PyFloat_FromDouble((double)timestep.gamma));
671
	    if(timestep.physical_time > -1)
Berk Onat's avatar
Berk Onat committed
672
	        PyDict_SetItemString(ret, "physical_time", PyFloat_FromDouble(timestep.physical_time));
673
674
675
676
677
678
	    //if(has_velocities > -1)
	    //    PyDict_SetItemString(ret, "has_velocities", PyLong_FromLong((long)has_velocities));
	    if(total_steps > -1)
	        PyDict_SetItemString(ret, "total_steps", PyLong_FromLong((long)total_steps));
	    PyDict_SetItemString(ret, "has_velocities", PyLong_FromLong((long)has_velocities));
            return (PyObject*) ret;
Berk Onat's avatar
Berk Onat committed
679
680
681
682
683
684
685
	}
    } else {
        PyErr_Format(PyExc_AttributeError, "molfile plugin does not have read_next_timestep function.");
	Py_RETURN_NONE;
    }
}

686
PyObject* are_plugins_same(PyObject* molpack_a, PyObject* molpack_b)
687
688
689
690
691
692
693
{
    molfile_plugin_t* plugin_a;
    molfile_plugin_t* plugin_b;
    PyObject *ret = NULL;
    MolObject* plugin_handle_a = (MolObject*) molpack_a;
    MolObject* plugin_handle_b = (MolObject*) molpack_b;
    if (plugin_handle_a->plugin) {
694
695
        plugin_a = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle_a->plugin);
        //plugin_a = plugin_handle_a->plugin;   
696
    } else {
697
        PyErr_Format(PyExc_IOError, "Arg 1 of the molfile plugin is not active.");
698
	return NULL;
699
700
    } 
    if (plugin_handle_b->plugin) {
701
702
        plugin_b = (molfile_plugin_t*) PyMolfileCapsule_AsVoidPtr(plugin_handle_b->plugin);
        //plugin_b = plugin_handle_b->plugin;   
703
    } else {
704
705
706
707
708
709
710
711
712
713
        PyErr_Format(PyExc_IOError, "Arg 2 of the molfile plugin is not active.");
	return NULL; 
    } 
    if(plugin_a == plugin_b){
      Py_RETURN_TRUE;
    } else {
      Py_RETURN_FALSE;
    }
}

714
PyObject* are_filehandles_same(PyObject* molpack_a, PyObject* molpack_b)
715
716
717
{
    MolObject* plugin_handle_a = (MolObject*) molpack_a;
    MolObject* plugin_handle_b = (MolObject*) molpack_b;
718
719
    void* file_handle_a; 
    void* file_handle_b; 
720
    if (plugin_handle_a->file_handle) {
721
722
        file_handle_a = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle_a->file_handle);
        //file_handle_a = plugin_handle_a->file_handle;   
723
    } else {
724
725
726
727
        PyErr_Format(PyExc_IOError, "no file handle in arg 1 of molfile plugin.");
	return NULL;
    } 
    if (plugin_handle_b->file_handle) {
728
729
        file_handle_b = (void*) PyMolfileCapsule_AsVoidPtr(plugin_handle_b->file_handle);
        //file_handle_b = plugin_handle_b->file_handle;   
730
731
    } else {
        PyErr_Format(PyExc_IOError, "no file handle in arg 2 of molfile plugin.");
732
733
734
735
736
737
738
739
740
	return NULL;
    } 
    if(file_handle_a == file_handle_b){
      Py_RETURN_TRUE;
    } else {
      Py_RETURN_FALSE;
    }
}