Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

py_atomsel.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 /***************************************************************************
00009  * RCS INFORMATION:
00010  *
00011  *      $RCSfile: py_atomsel.C,v $
00012  *      $Author: johns $        $Locker:  $             $State: Exp $
00013  *      $Revision: 1.49 $      $Date: 2020/07/23 03:27:52 $
00014  *
00015  ***************************************************************************
00016  * DESCRIPTION:
00017  *   New Python atom selection interface
00018  ***************************************************************************/
00019 #include "py_commands.h"
00020 #include "AtomSel.h"
00021 #include "VMDApp.h"
00022 #include "MoleculeList.h"
00023 #include "SymbolTable.h"
00024 #include "Measure.h"
00025 #include "SpatialSearch.h"
00026 
00027 // support for compiling against old (2.4.x and below) versions of Python
00028 #if PY_VERSION_HEX < ((2<<24)|(5<<16))
00029 typedef int Py_ssize_t;
00030 #endif
00031 
00032 typedef struct {
00033   PyObject_HEAD
00034   AtomSel *atomSel;
00035   VMDApp *app;
00036 } PyAtomSelObject;
00037 
00038 // Helper function to check if something is an atomsel type
00039 static int atomsel_Check(PyObject *obj) {
00040   if (PyObject_TypeCheck(obj, &Atomsel_Type))
00041     return 1;
00042   PyErr_SetString(PyExc_TypeError, "expected atomsel");
00043   return 0;
00044 }
00045 
00046 AtomSel *atomsel_AsAtomSel( PyObject *obj) {
00047   if (!atomsel_Check(obj))
00048     return NULL;
00049   return ((PyAtomSelObject *)obj)->atomSel;
00050 }
00051 
00052 
00053 // return molecule for atomsel, or NULL and set exception
00054 DrawMolecule *get_molecule(PyAtomSelObject *a) {
00055   int molid = a->atomSel->molid();
00056   DrawMolecule *mol = a->app->moleculeList->mol_from_id(molid);
00057   if (!mol) {
00058     PyErr_Format(PyExc_ValueError, "invalid molid '%d'", molid);
00059     return NULL;
00060   }
00061   return mol;
00062 }
00063 
00064 static const char atomsel_doc[] =
00065 "Create a new atom selection object\n\n"
00066 "Args:\n"
00067 "    selection (str): Atom selection string. Defaults to 'all'\n"
00068 "    molid (int): Molecule ID to select from. Defaults to -1 (top molecule)\n"
00069 "    frame (int): Frame to select on. Defaults to -1 (current)\n\n"
00070 "Example of usage:\n"
00071 "    >>> from vmd import atomsel\n"
00072 "    >>> s1 = atomsel('residue 1 to 10 and backbone')\n"
00073 "    >>> s1.resid\n"
00074 "     <snip> \n"
00075 "    >>> s1.beta = 5 # Set beta to 5 for all atoms in s1\n"
00076 "    >>> # Mass-weighted RMS alignment:\n"
00077 "    >>> mass = s1.get('mass')\n"
00078 "    >>> s2 = atomsel('residue 21 to 30 and backbone')\n"
00079 "    >>> mat = s1.fit(s2, mass)\n"
00080 "    >>> s1.move(mat)\n"
00081 "    >>> print s1.rmsd(s2)\n";
00082 
00083 // __del__(self)
00084 static void atomsel_dealloc( PyAtomSelObject *obj ) {
00085   delete obj->atomSel;
00086   ((PyObject *)(obj))->ob_type->tp_free((PyObject *)obj);
00087 }
00088 
00089 // __repr__(self)
00090 static PyObject *atomsel_repr(PyAtomSelObject *obj) {
00091   PyObject *result;
00092   AtomSel *sel;
00093   char *s;
00094 
00095   sel = obj->atomSel;
00096   s = new char[strlen(sel->cmdStr) + 100];
00097 
00098   sprintf(s, "atomsel('%s', molid=%d, frame=%d)", sel->cmdStr, sel->molid(),
00099           sel->which_frame);
00100   result = as_pystring(s);
00101 
00102   delete [] s;
00103   return result;
00104 }
00105 
00106 // __str__(self)
00107 static PyObject* atomsel_str(PyAtomSelObject *obj)
00108 {
00109   return as_pystring(obj->atomSel->cmdStr);
00110 }
00111 
00112 // __init__(self, selection, molid, frame)
00113 static PyObject *atomsel_new(PyTypeObject *type, PyObject *args,
00114                              PyObject *kwargs)
00115 {
00116   const char *kwlist[] = {"selection", "molid", "frame", NULL};
00117   int molid = -1, frame = -1;
00118   const char *sel = "all";
00119   DrawMolecule *mol;
00120   int parse_result;
00121   AtomSel *atomSel;
00122   PyObject *obj;
00123 
00124   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|sii:atomsel",
00125                                    (char**) kwlist, &sel, &molid, &frame))
00126     return NULL;
00127 
00128   VMDApp *app;
00129   if (!(app = get_vmdapp()))
00130     return NULL;
00131 
00132   if (molid < 0)
00133     molid = app->molecule_top();
00134 
00135   if (!valid_molid(molid, app))
00136     return NULL;
00137 
00138   mol = app->moleculeList->mol_from_id(molid);
00139 
00140   atomSel = new AtomSel(app, app->atomSelParser, mol->id());
00141   atomSel->which_frame = frame;
00142 
00143   Py_BEGIN_ALLOW_THREADS
00144     parse_result = atomSel->change(sel, mol);
00145   Py_END_ALLOW_THREADS
00146 
00147   if (parse_result == AtomSel::NO_PARSE) {
00148     PyErr_Format(PyExc_ValueError, "cannot parse atom selection text '%s'", sel);
00149     goto failure;
00150   }
00151 
00152   obj = type->tp_alloc(type, 0);
00153   if (!obj) {
00154     PyErr_Format(PyExc_MemoryError, "cannot allocate atomsel type");
00155     goto failure;
00156   }
00157 
00158   ((PyAtomSelObject *)obj)->app = app;
00159   ((PyAtomSelObject *)obj)->atomSel = atomSel;
00160 
00161   return obj;
00162 
00163 failure:
00164   delete atomSel;
00165   return NULL;
00166 }
00167 
00168 static const char get_doc[] =
00169 "Get attribute values for selected atoms\n\n"
00170 "Args:\n"
00171 "    attribute (str): Attribute to query\n"
00172 "Returns:\n"
00173 "    (list): Attribute value for each atom in selection";
00174 static PyObject *atomsel_get(PyAtomSelObject *a, PyObject *keyobj) {
00175   // FIXME: make everything in this pipeline const so that it's
00176   // thread-safe.
00177   const VMDApp *app = a->app;
00178   const AtomSel *atomSel = a->atomSel;
00179   const int num_atoms = a->atomSel->num_atoms;
00180 
00181   PyObject *newlist = NULL;
00182   SymbolTableElement *elem;
00183   SymbolTable *table;
00184   DrawMolecule *mol;
00185   int attrib_index, i;
00186   int *flgs;
00187   int j = 0;
00188 
00189   if (!(mol = get_molecule(a)))
00190     return NULL;
00191 
00192   const char *attr = as_constcharptr(keyobj);
00193   if (!attr) return NULL;
00194 
00195   //
00196   // Check for a valid attribute
00197   //
00198   table = app->atomSelParser;
00199   attrib_index = table->find_attribute(attr);
00200   if (attrib_index == -1) {
00201     PyErr_Format(PyExc_AttributeError, "unknown atom attribute '%s'", attr);
00202     return NULL;
00203   }
00204 
00205   elem = table->fctns.data(attrib_index);
00206   if (elem->is_a != SymbolTableElement::KEYWORD &&
00207       elem->is_a != SymbolTableElement::SINGLEWORD) {
00208     PyErr_Format(PyExc_AttributeError, "attribute '%s' is not a keyword or "
00209                  "singleword", attr);
00210     return NULL;
00211   }
00212 
00213   //
00214   // fetch the data
00215   //
00216   flgs = atomSel->on;
00217   atomsel_ctxt context(table, mol, atomSel->which_frame, attr);
00218 
00219   if (!(newlist = PyList_New(atomSel->selected)))
00220     return NULL;
00221 
00222   // Singleword attributes (backbone, etc) return array of bools
00223   if (elem->is_a == SymbolTableElement::SINGLEWORD) {
00224 
00225     int *tmp = new int[num_atoms];
00226     memcpy(tmp, atomSel->on, num_atoms*sizeof(int));
00227     elem->keyword_single(&context, num_atoms, tmp);
00228 
00229     for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00230       if (flgs[i]) {
00231         PyObject *val = tmp[i] ? Py_True : Py_False;
00232         Py_INCREF(val); // SET_ITEM steals a reference so we make one to steal
00233         PyList_SET_ITEM(newlist, j++, val);
00234 
00235         if (PyErr_Occurred())
00236           goto failure;
00237       }
00238     }
00239     delete [] tmp;
00240 
00241   // String attributes
00242   } else if (table->fctns.data(attrib_index)->returns_a \
00243           == SymbolTableElement::IS_STRING) {
00244 
00245     const char **tmp= new const char *[num_atoms];
00246     elem->keyword_string(&context, num_atoms, tmp, flgs);
00247 
00248     for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00249       if (flgs[i]) {
00250         PyList_SET_ITEM(newlist, j++, as_pystring(tmp[i]));
00251 
00252         if (PyErr_Occurred()) {
00253           delete [] tmp;
00254           goto failure;
00255         }
00256       }
00257     }
00258     delete [] tmp;
00259 
00260   // Integer attributes
00261   } else if (table->fctns.data(attrib_index)->returns_a \
00262           == SymbolTableElement::IS_INT) {
00263 
00264     int *tmp = new int[num_atoms];
00265     elem->keyword_int(&context, num_atoms, tmp, flgs);
00266 
00267     for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00268       if (flgs[i]) {
00269         PyList_SET_ITEM(newlist, j++, as_pyint(tmp[i]));
00270 
00271         if (PyErr_Occurred()) {
00272           delete [] tmp;
00273           goto failure;
00274         }
00275       }
00276     }
00277     delete [] tmp;
00278 
00279   // Floating point attributes
00280   } else if (table->fctns.data(attrib_index)->returns_a \
00281           == SymbolTableElement::IS_FLOAT) {
00282 
00283     double *tmp = new double[num_atoms];
00284     elem->keyword_double(&context, num_atoms, tmp, flgs);
00285 
00286     for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00287       if (flgs[i]) {
00288         PyList_SET_ITEM(newlist, j++, PyFloat_FromDouble(tmp[i]));
00289 
00290         if (PyErr_Occurred()) {
00291           delete [] tmp;
00292           goto failure;
00293         }
00294 
00295       }
00296     }
00297     delete [] tmp;
00298   }
00299 
00300   return newlist;
00301 
00302 failure:
00303   return NULL;
00304 }
00305 
00306 static const char listattrs_doc[] =
00307 "List available atom attributes\n"
00308 "Args:\n"
00309 "    changeable (bool): If only user-changeable attributes should be listed\n"
00310 "        Defaults to False\n"
00311 "Returns:\n"
00312 "    (list of str): Atom attributes. These attributes may be accessed or\n"
00313 "        set with a . after the class name.\n"
00314 "Example to list available attributes and get the x coordinate attribute:\n"
00315 "    >>> sel = atomsel('protein')\n"
00316 "    >>> sel.list_attributes()\n"
00317 "    >>> sel.x";
00318 static PyObject *py_list_attrs(PyAtomSelObject *obj, PyObject *args,
00319                                PyObject *kwargs)
00320 {
00321   const char *kwlist[] = {"changeable", NULL};
00322   SymbolTableElement *elem;
00323   PyObject *result = NULL;
00324   int changeable = 0;
00325   SymbolTable *table;
00326   int i;
00327 
00328   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|O&:atomsel.list_attributes",
00329                                    (char**) kwlist, convert_bool, &changeable))
00330     return NULL;
00331 
00332   table = obj->app->atomSelParser;
00333   result = PyList_New(0);
00334 
00335   if (!result)
00336     goto failure;
00337 
00338   for (i = 0; i< table->fctns.num(); i++) {
00339     // Only list singleword or keyword attributes that are gettable
00340     elem = table->fctns.data(i);
00341     if (elem->is_a != SymbolTableElement::KEYWORD
00342      && elem->is_a != SymbolTableElement::SINGLEWORD)
00343       continue;
00344 
00345     // Only list changeable attributes if requested
00346     if (changeable && !table->is_changeable(i))
00347       continue;
00348 
00349     // Don't list unnamed attributes
00350     if (!table->fctns.name(i) || !strlen(table->fctns.name(i)))
00351       continue;
00352 
00353     // Don't leak references with PyList_Append
00354     PyObject *attr = as_pystring(table->fctns.name(i));
00355     PyList_Append(result, attr);
00356     Py_XDECREF(attr);
00357 
00358     if (PyErr_Occurred())
00359       goto failure;
00360   }
00361 
00362   return result;
00363 
00364 failure:
00365   PyErr_SetString(PyExc_RuntimeError, "Problem listing attributes");
00366   Py_XDECREF(result);
00367   return NULL;
00368 }
00369 
00370 static PyObject *legacy_atomsel_get(PyObject *o, PyObject *args,
00371                                     PyObject *kwargs)
00372 {
00373   const char *kwlist[] = {"attribute", NULL};
00374   PyObject *attr, *result;
00375 
00376   PyErr_Warn(PyExc_DeprecationWarning, "atomsel.get is deprecated. You can\n"
00377              "directly query the attributes of this object instead");
00378 
00379   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O:atomsel.get",
00380                                    (char**) kwlist, &attr))
00381     return NULL;
00382 
00383   // It sets an exception on error
00384   if (!(result = atomsel_get((PyAtomSelObject*) o, attr)))
00385     return NULL;
00386 
00387   return result;
00388 }
00389 
00390 // __getattr__, with support for defined methods too
00391 static PyObject *atomsel_getattro(PyObject *o, PyObject *attr_name)
00392 {
00393   // Check that there isn't a class method by this name first
00394   PyObject *tmp;
00395   if (!(tmp = PyObject_GenericGetAttr(o, attr_name))) {
00396     // If no method found, throws AttributeError that we clear.
00397     // If it threw a different exception though, we do fail.
00398     if (!PyErr_ExceptionMatches(PyExc_AttributeError))
00399       return NULL;
00400     PyErr_Clear();
00401   } else {
00402     return tmp;
00403   }
00404 
00405   return atomsel_get((PyAtomSelObject*) o, attr_name);
00406 }
00407 
00408 static void help_int(void* p, PyObject *obj)
00409 {
00410   *((int*)p) = as_int(obj);
00411 }
00412 
00413 static void help_double(void* p, PyObject *obj)
00414 {
00415   *((double*)p) = PyFloat_AsDouble(obj);
00416 }
00417 
00418 static void help_constcharptr(void* p, PyObject *obj)
00419 {
00420   const char *str = as_constcharptr(obj);
00421   *((const char**)p) = str;
00422 }
00423 
00424 // Helper functions for set
00425 static int build_set_values(const void* list, const AtomSel *atomSel, 
00426                             PyObject* val, int dtype_size,
00427                             void (*converter)(void*, PyObject*)) {
00428   int i, is_array;
00429   int j = 0;
00430   void *p;
00431 
00432   // Determine if passed PyObject is an array
00433   is_array = (PySequence_Check(val) && !is_pystring(val)) ? 1 : 0;
00434 
00435   const int *flgs = atomSel->on;
00436   for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00437     // Continue if atom is not part of selection
00438     if (!flgs[i])
00439       continue;
00440 
00441     // Get the correct atom
00442     p = (char*) list + i*dtype_size;
00443 
00444     // Set ival to each value in array if we have an array
00445     // If array has length 1, repeat first array element*
00446     // *I think this shouldn't be supported, but keeping it for compatilibity
00447     if (is_array) {
00448       if (PySequence_Length(val) == 1)
00449         converter(p, PySequence_ITEM(val, 0));
00450       else
00451         converter(p, PySequence_ITEM(val, j++));
00452 
00453     // If not an array of values, unpack PyObject alone
00454     } else {
00455       converter(p, val);
00456     }
00457     if (PyErr_Occurred())
00458       goto failure;
00459   }
00460 
00461   return 0;
00462 
00463 failure:
00464   PyErr_SetString(PyExc_ValueError, "sequence passed to set contains "
00465                   "the wrong type of data (integer, float, or str)");
00466   return 1;
00467 }
00468 
00469 
00470 static int atomsel_set(PyAtomSelObject *a, const char *name, PyObject *val)
00471 {
00472   const VMDApp *app = a->app;
00473   const AtomSel *atomSel = a->atomSel;
00474   const int num_atoms = atomSel->num_atoms;
00475   SymbolTable *table = app->atomSelParser;
00476   int *flgs = atomSel->on;
00477   void* list = NULL;
00478 
00479   SymbolTableElement *elem;
00480   DrawMolecule *mol;
00481   int attrib_index;
00482 
00483   // Check for a valid molecule
00484   if (!(mol = get_molecule(a))) {
00485     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
00486     return -1;
00487   }
00488 
00489   // Check for a valid attribute
00490   if ((attrib_index = table->find_attribute(name)) == -1) {
00491     PyErr_Format(PyExc_AttributeError, "unknown atom attribute '%s'", name);
00492     return -1;
00493   }
00494 
00495   elem = table->fctns.data(attrib_index);
00496   if (elem->is_a != SymbolTableElement::KEYWORD &&
00497       elem->is_a != SymbolTableElement::SINGLEWORD) {
00498     PyErr_Format(PyExc_AttributeError, "attribute '%s' is not a keyword or "
00499                  "singleword", name);
00500     return -1;
00501   }
00502 
00503   if (!table->is_changeable(attrib_index)) {
00504     PyErr_Format(PyExc_AttributeError, "attribute '%s' is not modifiable",
00505                  name);
00506     return -1;
00507   }
00508 
00509   // singlewords can never be set, so macro is NULL.
00510   atomsel_ctxt context(table, mol, atomSel->which_frame, NULL);
00511 
00512   // Integer type
00513   if (elem->returns_a == SymbolTableElement::IS_INT) {
00514     list = malloc(num_atoms * sizeof(int));
00515     if (build_set_values(list, atomSel, val, sizeof(int), help_int))
00516       goto failure;
00517 
00518     elem->set_keyword_int(&context, num_atoms, (int*) list, atomSel->on);
00519 
00520   // Float type
00521   } else if (elem->returns_a == SymbolTableElement::IS_FLOAT) {
00522     list = malloc(num_atoms * sizeof(double));
00523     if (build_set_values(list, atomSel, val, sizeof(double), help_double))
00524       goto failure;
00525 
00526     elem->set_keyword_double(&context, num_atoms, (double*) list, flgs);
00527 
00528   } else if (elem->returns_a == SymbolTableElement::IS_STRING) {
00529     list = malloc(num_atoms * sizeof(char*));
00530     if (build_set_values(list, atomSel, val, sizeof(char*), help_constcharptr))
00531         goto failure;
00532 
00533     elem->set_keyword_string(&context, num_atoms, (const char**) list, atomSel->on);
00534     // No special memory management needed for strings as python API
00535     // as_constcharptr gives a pointer into PyObject* that Python is managing
00536   }
00537 
00538   free(list);
00539   list = NULL;
00540 
00541   // Recompute the color assignments if certain atom attributes are changed.
00542   if (!strcmp(name, "name") ||
00543       !strcmp(name, "type") ||
00544       !strcmp(name, "resname") ||
00545       !strcmp(name, "chain") ||
00546       !strcmp(name, "segid") ||
00547       !strcmp(name, "segname"))
00548     app->moleculeList->add_color_names(mol->id());
00549 
00550   mol->force_recalc(DrawMolItem::SEL_REGEN | DrawMolItem::COL_REGEN);
00551 
00552   return 0;
00553 
00554 failure:
00555   // Exception already set by build_set_values or similar
00556   free(list);
00557   return 1;
00558 }
00559 
00560 // Legacy atomsel.set method
00561 static const char set_doc[] =
00562 "Set attributes for selected atoms\n\n"
00563 "Args:\n"
00564 "    attribute (str): Attribute to set\n"
00565 "    value: Value for attribute. If single value, all atoms will be set to\n"
00566 "        have the same value. Otherwise, pass a sequence (list or tuple) of\n"
00567 "        values, one per atom in selection";
00568 static PyObject* legacy_atomsel_set(PyObject *o, PyObject *args,
00569                                     PyObject *kwargs)
00570 {
00571   const char *kwlist[] = {"attribute", "value", NULL};
00572   const char *attr;
00573   PyObject *val;
00574 
00575   PyErr_Warn(PyExc_DeprecationWarning, "atomsel.set is deprecated. You can\n"
00576              "directly assign to the attributes instead.");
00577 
00578   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "sO:atomsel.set",
00579                                   (char**) kwlist,  &attr, &val))
00580     return NULL;
00581 
00582   // It sets an exception on error
00583   if (atomsel_set((PyAtomSelObject*) o, attr, val))
00584     return NULL;
00585 
00586   Py_INCREF(Py_None);
00587   return Py_None;
00588 }
00589 
00590 // __setattr__, with support for defined methods too
00591 static int atomsel_setattro(PyObject *o, PyObject *name, PyObject *value)
00592 {
00593   // Check that there isn't a class method for this set first
00594   // GenericSetAttr returns 0 on success
00595   // If no setter method found, throws AttributeError that we clear
00596   // We do fail if a different exception is thrown
00597   if (PyObject_GenericSetAttr(o, name, value)) {
00598     if (!PyErr_ExceptionMatches(PyExc_AttributeError))
00599       return -1;
00600     PyErr_Clear();
00601 
00602   } else {
00603     return 0;
00604   }
00605 
00606   return atomsel_set((PyAtomSelObject*) o, as_constcharptr(name), value);
00607 }
00608 
00609 static const char frame_doc[] =
00610 "Get the frame an atomsel object references. Changing the frame does not\n"
00611 "immediately update the selection. Use `atomsel.update()` to do that.\n"
00612 "Special frame values are -1 for the current frame and -2 for the last frame";
00613 static PyObject *getframe(PyAtomSelObject *a, void *) {
00614 
00615   AtomSel *atomSel = a->atomSel;
00616   DrawMolecule *mol;
00617 
00618   if (!(mol = get_molecule(a)))
00619     return NULL;
00620 
00621   return as_pyint(atomSel->which_frame);
00622 }
00623 
00624 static int setframe(PyAtomSelObject *a, PyObject *frameobj, void *) {
00625 
00626   AtomSel *atomSel = a->atomSel;
00627   DrawMolecule *mol;
00628   int frame;
00629   if (!(mol = get_molecule(a)))
00630     return -1;
00631 
00632   frame = as_int(frameobj);
00633   if (PyErr_Occurred())
00634     return -1;
00635 
00636   if (frame != AtomSel::TS_LAST && frame != AtomSel::TS_NOW &&
00637       (frame < 0 || frame >= mol->numframes())) {
00638     PyErr_Format(PyExc_ValueError, "Frame '%d' invalid for this molecule",
00639                  frame);
00640     return -1;
00641   }
00642 
00643   atomSel->which_frame = frame;
00644   return 0;
00645 }
00646 
00647 static const char update_doc[] =
00648 "Recompute which atoms in the molecule belong to this selection. For example\n"
00649 "when the selection is distance base (e.g. it uses 'within'), changes to the\n"
00650 "frame of this atom selection will not be reflected in the selected atoms\n"
00651 "until this method is called";
00652 static PyObject *py_update(PyAtomSelObject *a) {
00653 
00654   AtomSel *atomSel = a->atomSel;
00655   DrawMolecule *mol;
00656   if (!(mol = get_molecule(a))) return NULL;
00657 
00658   Py_BEGIN_ALLOW_THREADS
00659     atomSel->change(NULL, mol);
00660   Py_END_ALLOW_THREADS
00661 
00662   Py_INCREF(Py_None);
00663   return Py_None;
00664 }
00665 
00666 static const char write_doc[] =
00667 "Write atoms in this selection to a file of the given type\n\n"
00668 "Args:\n"
00669 "    filetype (str): File type to write, as defined by molfileplugin\n"
00670 "    filename (str): Filename to write to";
00671 static PyObject *py_write(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
00672 {
00673   const char *kwlist[] = {"filetype", "filename", NULL};
00674   const char *filetype, *filename;
00675   AtomSel *atomSel = a->atomSel;
00676   DrawMolecule *mol;
00677   int frame = -1;
00678   FileSpec spec;
00679   VMDApp *app;
00680   int rc;
00681 
00682   if (!(mol = get_molecule(a))) {
00683     PyErr_SetString(PyExc_ValueError, "molecule for this selection is deleted");
00684     return NULL;
00685   }
00686 
00687   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ss:atomsel.write",
00688                                    (char**) kwlist, &filetype, &filename))
00689     return NULL;
00690 
00691   if (!(app = get_vmdapp()))
00692     return NULL;
00693 
00694   switch (atomSel->which_frame) {
00695     case AtomSel::TS_NOW:  frame = mol->frame(); break;
00696     case AtomSel::TS_LAST: frame = mol->numframes()-1; break;
00697     default:               frame = atomSel->which_frame; break;
00698   }
00699 
00700   // If frame out of bounds, return error
00701   // (formerly this just truncated frame to being in the valid range)
00702   if (frame < 0 || frame >= mol->numframes()) {
00703     PyErr_Format(PyExc_ValueError, "frame '%d' out of bounds for this molecule",
00704                  frame);
00705     return NULL;
00706   }
00707 
00708   // Write the requested atoms to the file
00709   spec.first = frame;
00710   spec.last = frame;
00711   spec.stride = 1;
00712   spec.waitfor = FileSpec::WAIT_ALL;
00713   spec.selection = atomSel->on;
00714 
00715   Py_BEGIN_ALLOW_THREADS
00716     rc = app->molecule_savetrajectory(mol->id(), filename, filetype, &spec);
00717   Py_END_ALLOW_THREADS
00718 
00719   if (rc < 0) {
00720     PyErr_SetString(PyExc_ValueError, "Unable to write selection to file.");
00721     return NULL;
00722   }
00723 
00724   Py_INCREF(Py_None);
00725   return Py_None;
00726 }
00727 
00728 static const char bonds_doc[] =
00729 "For each atom in selection, a list of the indices of atoms to which it is\n"
00730 "bonded.\nTo set bonds, pass a sequence with length of the number of atoms in\n"
00731 "the selection, with each entry in the sequence being a list or tuple\n"
00732 "containing the atom indices to which that atom in the selection is bound\n\n"
00733 "For example, for a water molecule with atoms H-O-H:\n"
00734 ">>> sel = atomsel('water and residue 0')\n"
00735 ">>> sel.bonds = [(1), (0,2), (1)]";
00736 static PyObject *getbonds(PyAtomSelObject *a, void *)
00737 {
00738   PyObject *bondlist = NULL, *newlist = NULL;
00739   AtomSel *atomSel = a->atomSel;
00740   DrawMolecule *mol;
00741   int i, j, k = 0;
00742 
00743   if (!(mol = get_molecule(a))) {
00744     PyErr_SetString(PyExc_ValueError, "selection is on deleted molecule");
00745     return NULL;
00746   }
00747 
00748   newlist = PyList_New(atomSel->selected);
00749   if (!newlist)
00750     goto failure;
00751 
00752   for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00753 
00754     if (!atomSel->on[i])
00755       continue;
00756 
00757     const MolAtom *atom = mol->atom(i);
00758     bondlist = PyList_New(atom->bonds);
00759     if (!bondlist)
00760       goto failure;
00761 
00762     for (j=0; j<atom->bonds; j++) {
00763       PyList_SET_ITEM(bondlist, j, as_pyint(atom->bondTo[j]));
00764       if (PyErr_Occurred())
00765         goto failure;
00766     }
00767 
00768     PyList_SET_ITEM(newlist, k++, bondlist);
00769     if (PyErr_Occurred())
00770       goto failure;
00771   }
00772   return newlist;
00773 
00774 failure:
00775   Py_XDECREF(newlist);
00776   Py_XDECREF(bondlist);
00777   PyErr_SetString(PyExc_RuntimeError, "Problem getting bonds");
00778   return NULL;
00779 }
00780 
00781 static int setbonds(PyAtomSelObject *a, PyObject *obj, void *)
00782 {
00783   AtomSel *atomSel = a->atomSel;
00784   DrawMolecule *mol;
00785   PyObject *atomids;
00786   int ibond = 0;
00787   int i, j, k;
00788 
00789   if (!(mol = get_molecule(a))) {
00790     PyErr_SetString(PyExc_ValueError, "selection is on deleted molecule");
00791     return -1;
00792   }
00793 
00794   if (!PySequence_Check(obj)) {
00795     PyErr_SetString(PyExc_TypeError, "Argument to setbonds must be a sequence");
00796     return -1;
00797   }
00798 
00799   if (PySequence_Size(obj) != atomSel->selected) {
00800     PyErr_SetString(PyExc_ValueError, "atomlist and bondlist must be the same "
00801                     "size");
00802     return -1;
00803   }
00804 
00805   mol->force_recalc(DrawMolItem::MOL_REGEN); // many reps ignore bonds
00806   for (i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
00807     if (!atomSel->on[i])
00808       continue;
00809 
00810     MolAtom *atom = mol->atom(i);
00811     if (!(atomids = PySequence_GetItem(obj, ibond++))) {
00812       PyErr_Format(PyExc_RuntimeError, "Could not get bonds for atom %d",
00813                    ibond - 1);
00814       goto failure;
00815     }
00816 
00817     if (!PySequence_Check(atomids)) {
00818       PyErr_SetString(PyExc_TypeError, "Bonded atoms must be a sequence");
00819       return -1;
00820     }
00821 
00822     k = 0;
00823     for (j = 0; j < PySequence_Size(atomids); j++) {
00824       int bond = as_int(PySequence_GetItem(atomids, j));
00825 
00826       if (PyErr_Occurred()) {
00827         PyErr_SetString(PyExc_TypeError, "atom indices to set bonds must be "
00828                         "integers");
00829         goto failure;
00830       }
00831 
00832       if (bond >= 0 && bond < mol->nAtoms) {
00833         atom->bondTo[k++] = bond;
00834       }
00835     }
00836     atom->bonds = k;
00837   }
00838 
00839   return 0;
00840 
00841 failure:
00842   return -1;
00843 }
00844 
00845 static const char molid_doc[] =
00846 "The molecule ID the selection is associated with";
00847 static PyObject *getmolid(PyAtomSelObject *a, void *) {
00848   if (!a->app->moleculeList->mol_from_id(a->atomSel->molid())) {
00849     PyErr_SetString(PyExc_ValueError, "selection is on deleted molecule");
00850     return NULL;
00851   }
00852   return as_pyint(a->atomSel->molid());
00853 }
00854 
00855 static PyGetSetDef atomsel_getset[] = {
00856   {(char*) "frame", (getter)getframe, (setter)setframe, (char*) frame_doc, NULL},
00857   {(char*) "bonds", (getter)getbonds, (setter)setbonds, (char*) bonds_doc, NULL},
00858   {(char*) "molid", (getter)getmolid, (setter)NULL, (char*) molid_doc, NULL},
00859   {NULL },
00860 };
00861 
00862 // utility routine for parsing weight values.  Uses the sequence protocol
00863 // so that sequence-type structure (list, tuple) will be accepted.
00864 static float *parse_weight(AtomSel *sel, PyObject *wtobj)
00865 {
00866   float *weight = new float[sel->selected];
00867   PyObject *seq = NULL;
00868   int i;
00869 
00870   // If no weights passed, set them all to 1.0
00871   if (!wtobj || wtobj == Py_None) {
00872     for (int i=0; i<sel->selected; i++)
00873       weight[i] = 1.0f;
00874     return weight;
00875   }
00876 
00877   if (!(seq = PySequence_Fast(wtobj, "weight must be a sequence.")))
00878     goto failure;
00879 
00880   if (PySequence_Size(seq) != sel->selected) {
00881     PyErr_SetString(PyExc_ValueError, "weight must be same size as selection.");
00882     goto failure;
00883   }
00884 
00885   for (i = 0; i < sel->selected; i++) {
00886     PyObject *w = PySequence_Fast_GET_ITEM(seq, i);
00887     if (!PyFloat_Check(w)) {
00888       PyErr_SetString(PyExc_TypeError, "Weights must be floating point numbers");
00889       goto failure;
00890     }
00891 
00892     double tmp = PyFloat_AsDouble(PySequence_Fast_GET_ITEM(seq, i));
00893     weight[i] = (float)tmp;
00894   }
00895 
00896   Py_XDECREF(seq);
00897   return weight;
00898 
00899 failure:
00900   Py_XDECREF(seq);
00901   delete [] weight;
00902   return NULL;
00903 }
00904 
00905 static const char minmax_doc[] =
00906 "Get minimum and maximum coordinates for selected atoms\n\n"
00907 "Args:\n"
00908 "    radii (bool): If atomic radii should be included in calculation\n"
00909 "        Defaults to False.\n"
00910 "Returns:\n"
00911 "    (2-tuple of tuples): (x,y,z) coordinate of minimum, then maximum atom";
00912 static PyObject *minmax(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
00913 {
00914   const char *kwlist[] = {"radii", NULL};
00915   AtomSel *atomSel = a->atomSel;
00916   const float *radii;
00917   float min[3], max[3];
00918   DrawMolecule *mol;
00919   int withradii = 0;
00920   const float *pos;
00921   int rc;
00922 
00923   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|O&:atomsel.minmax",
00924                                    (char**) kwlist, convert_bool, &withradii))
00925     return NULL;
00926 
00927   if (!(mol = get_molecule(a))) {
00928     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
00929     return NULL;
00930   }
00931 
00932   radii = withradii ? mol->extraflt.data("radius") : NULL;
00933 
00934   pos = atomSel->coordinates(a->app->moleculeList);
00935   rc = measure_minmax(atomSel->num_atoms, atomSel->on, pos, radii, min, max);
00936   if (rc < 0) {
00937     PyErr_SetString(PyExc_ValueError, measure_error(rc));
00938     return NULL;
00939   }
00940 
00941   return Py_BuildValue("(f,f,f),(f,f,f)",
00942       min[0], min[1], min[2], max[0], max[1], max[2]);
00943 }
00944 
00945 static const char centerperres_doc[] =
00946 "Get the coordinates of the center of each residue in the selection,\n"
00947 "optionally weighted by weight\n\n"
00948 "Args:\n"
00949 "    weight (list of float): Weights for each atom in selection. Optional.\n"
00950 "        weights cannot be 0 otherwise NaN will be returned.\n"
00951 "Returns:\n"
00952 "    (list of tuple): (x,y,z) coordinates of center of each residue";
00953 static PyObject *centerperresidue(PyAtomSelObject *a, PyObject *args,
00954                                   PyObject *kwargs)
00955 {
00956   PyObject *weightobj = NULL, *returnlist = NULL, *element = NULL;
00957   const char *kwlist[] = {"weight", NULL};
00958   AtomSel *atomSel = a->atomSel;
00959   float *weight, *cen;
00960   DrawMolecule *mol;
00961   int i, j, ret_val;
00962 
00963   if (!(mol = get_molecule(a))) {
00964     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
00965     return NULL;
00966   }
00967 
00968   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|O:atomsel.centerperresidue",
00969                                    (char**) kwlist, &weightobj))
00970     return NULL;
00971 
00972   weight = parse_weight(atomSel, weightobj);
00973   if (!weight)
00974     return NULL;
00975 
00976   // compute center
00977   cen = new float[3*atomSel->selected];
00978   ret_val = measure_center_perresidue(a->app->moleculeList, atomSel,
00979                                       atomSel->coordinates(a->app->moleculeList),
00980                                       weight, cen);
00981   delete [] weight;
00982 
00983   if (ret_val < 0) {
00984     PyErr_SetString(PyExc_ValueError, measure_error(ret_val));
00985     goto failure;
00986     return NULL;
00987   }
00988 
00989   //Build the python list.
00990   returnlist = PyList_New(ret_val);
00991   for (i = 0; i < ret_val; i++) {
00992     element = PyTuple_New(3);
00993 
00994     for (j = 0; j < 3; j++) {
00995       PyTuple_SET_ITEM(element, j, PyFloat_FromDouble((double) cen[3*i+j]));
00996       if (PyErr_Occurred()) {
00997         PyErr_SetString(PyExc_ValueError, "Problem building center list");
00998         goto failure;
00999       }
01000     }
01001 
01002     PyList_SET_ITEM(returnlist, i, element);
01003     if (PyErr_Occurred()) {
01004       PyErr_SetString(PyExc_ValueError, "Problem building center list");
01005       goto failure;
01006     }
01007   }
01008 
01009   delete [] cen;
01010   return returnlist;
01011 
01012 failure:
01013   delete [] cen;
01014   delete [] weight;
01015   Py_XDECREF(returnlist);
01016   Py_XDECREF(element);
01017   return NULL;
01018 }
01019 
01020 
01021 static const char rmsfperres_doc[] =
01022 "Measures the root-mean-square fluctuation (RMSF) along a trajectory on a\n"
01023 "per-residue basis. RMSF is the mean deviation from the average position\n\n"
01024 "Args:\n"
01025 "    first (int): First frame to include. Defaults to 0 (beginning).\n"
01026 "    last (int): Last frame to include. Defaults to -1 (end).\n"
01027 "    step (int): Use every step-th frame. Defaults to 1 (all frames)\n"
01028 "Returns:\n"
01029 "    (list of float): RMSF for each residue in selection";
01030 static PyObject *py_rmsfperresidue(PyAtomSelObject *a, PyObject *args,
01031                                    PyObject *kwargs)
01032 {
01033   const char *kwlist[] = {"first", "last", "step", NULL};
01034   int first=0, last=-1, step=1;
01035   PyObject *returnlist = NULL;
01036   int ret_val, i;
01037   float *rmsf;
01038 
01039   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|iii:atomsel.rmsfperresidue",
01040                                    (char**) kwlist, &first, &last, &step))
01041     return NULL;
01042 
01043   // Check molecule is valid
01044   if (!get_molecule(a)) {
01045     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01046     return NULL;
01047   }
01048 
01049   rmsf = new float[a->atomSel->selected];
01050   ret_val = measure_rmsf_perresidue(a->atomSel, a->app->moleculeList,
01051                                     first, last, step, rmsf);
01052   if (ret_val < 0) {
01053     PyErr_SetString(PyExc_RuntimeError, measure_error(ret_val));
01054     goto failure;
01055   }
01056 
01057   //Build the python list.
01058   returnlist = PyList_New(ret_val);
01059   for (i = 0; i < ret_val; i++) {
01060     PyList_SetItem(returnlist, i, PyFloat_FromDouble(rmsf[i]));
01061     if (PyErr_Occurred()) {
01062       PyErr_SetString(PyExc_RuntimeError, "Problem building rmsf list");
01063       goto failure;
01064     }
01065   }
01066 
01067   delete [] rmsf;
01068   return returnlist;
01069 
01070 failure:
01071   Py_XDECREF(returnlist);
01072   delete [] rmsf;
01073   return NULL;
01074 }
01075 
01076 static const char center_doc[] =
01077 "Get the coordinates of the center of the selection, optionally with weights\n"
01078 "on the selection atoms\n\n"
01079 "Args:\n"
01080 "    weight (list of float): Weight on each atom. Optional\n"
01081 "Returns:\n"
01082 "    (3-tuple of float): (x,y,z) coordinates of center of the selection";
01083 static PyObject *center(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01084 {
01085   const char *kwlist[] = {"weight", NULL};
01086   AtomSel *atomSel = a->atomSel;
01087   PyObject *weightobj = NULL;
01088   DrawMolecule *mol;
01089   float *weight;
01090   float cen[3];
01091   int ret_val;
01092 
01093   if (!(mol = get_molecule(a))) {
01094     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01095     return NULL;
01096   }
01097 
01098   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|O:atomsel.center",
01099                                    (char**) kwlist, &weightobj))
01100     return NULL;
01101 
01102   weight = parse_weight(atomSel, weightobj);
01103   if (!weight)
01104     return NULL;
01105 
01106   // compute center
01107   ret_val = measure_center(atomSel, atomSel->coordinates(a->app->moleculeList),
01108                            weight, cen);
01109   delete [] weight;
01110 
01111   if (ret_val < 0) {
01112     PyErr_SetString(PyExc_ValueError, measure_error(ret_val));
01113     return NULL;
01114   }
01115 
01116   return Py_BuildValue("(f,f,f)", cen[0], cen[1], cen[2]);
01117 }
01118 
01119 // helper function to validate weights with multiple atomsels
01120 static float *parse_two_selections_return_weight(PyAtomSelObject *a,
01121                                                  PyObject *args,
01122                                                  PyObject *kwargs,
01123                                                  AtomSel **othersel)
01124 {
01125   const char *kwlist[] = {"selection", "weight", NULL};
01126   AtomSel *atomSel = a->atomSel;
01127   AtomSel *sel2;
01128   PyObject *weightobj = NULL;
01129   PyObject *other;
01130   DrawMolecule *mol;
01131   float *weight;
01132 
01133   if (!(mol = get_molecule(a))){
01134     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01135     return NULL;
01136   }
01137 
01138   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!|O", (char**) kwlist,
01139                                    &Atomsel_Type, &other, &weightobj))
01140     return NULL;
01141 
01142   weight = parse_weight(atomSel, weightobj);
01143   if (!weight)
01144     return NULL;
01145 
01146   sel2 = ((PyAtomSelObject *)other)->atomSel;
01147   if (!get_molecule((PyAtomSelObject *)other)) {
01148     PyErr_SetString(PyExc_ValueError,
01149                     "selection argument is on a deleted molecule");
01150     delete [] weight;
01151     return NULL;
01152   }
01153 
01154   if (atomSel->selected != sel2->selected) {
01155     PyErr_SetString(PyExc_ValueError, "Selections must have same number of "
01156                     "atoms");
01157     delete [] weight;
01158     return NULL;
01159   }
01160 
01161   *othersel = sel2;
01162   return weight;
01163 }
01164 
01165 static const char rmsdperres_doc[] =
01166 "Get the per-residue root-mean-square (RMS) distance between selections\n\n"
01167 "Args:\n"
01168 "    selection (atomsel): Selection to compute distance to. Must have the\n"
01169 "        same number of atoms as this selection\n"
01170 "    weight (list of float): Weight for atoms, one per atom in selection.\n"
01171 "        Optional\n"
01172 "Returns:\n"
01173 "    (list of float): RMSD between each residue in selection";
01174 static PyObject *py_rmsdperresidue(PyAtomSelObject *a, PyObject *args,
01175                                    PyObject *kwargs)
01176 {
01177   PyObject *returnlist = NULL;
01178   float *weight, *rmsd;
01179   AtomSel *sel2;
01180   int i, rc;
01181 
01182   weight = parse_two_selections_return_weight(a, args, kwargs, &sel2);
01183   if (!weight)
01184     return NULL;
01185 
01186   rmsd = new float[a->atomSel->selected];
01187   rc = measure_rmsd_perresidue(a->atomSel, sel2, a->app->moleculeList,
01188                                a->atomSel->selected, weight, rmsd);
01189   delete [] weight;
01190 
01191   if (rc < 0) {
01192     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01193     goto failure;
01194   }
01195 
01196   //Build the python list.
01197   returnlist = PyList_New(rc);
01198   for (i = 0; i < rc; i++) {
01199     PyList_SET_ITEM(returnlist, i, PyFloat_FromDouble(rmsd[i]));
01200     if (PyErr_Occurred()) {
01201       PyErr_SetString(PyExc_RuntimeError, "Problem building rmsd list");
01202       goto failure;
01203     }
01204   }
01205 
01206   delete [] rmsd;
01207   return returnlist;
01208 
01209 failure:
01210   delete [] rmsd;
01211   Py_XDECREF(returnlist);
01212   return NULL;
01213 }
01214 
01215 static const char rmsd_doc[] =
01216 "Calculate the root-mean-square distance (RMSD) between selections. Atoms\n"
01217 "must be in the same order in each selection\n\n"
01218 "Args:\n"
01219 "    selection (atomsel): Other selection to compute RMSD to. Must have\n"
01220 "        the same number of atoms as this selection\n"
01221 "    weight (list of float): Weight per atom, optional\n"
01222 "Returns:\n"
01223 "    (float): RMSD between selections.";
01224 static PyObject *py_rmsd(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01225 {
01226   AtomSel *sel2;
01227   float *weight;
01228   float rmsd;
01229   int rc;
01230 
01231   weight = parse_two_selections_return_weight(a, args, kwargs, &sel2);
01232   if (!weight)
01233     return NULL;
01234 
01235   rc = measure_rmsd(a->atomSel, sel2, a->atomSel->selected,
01236                     a->atomSel->coordinates(a->app->moleculeList),
01237                     sel2->coordinates(a->app->moleculeList),
01238                     weight, &rmsd);
01239   delete [] weight;
01240 
01241   if (rc < 0) {
01242     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01243     return NULL;
01244   }
01245 
01246   return PyFloat_FromDouble(rmsd);
01247 }
01248 
01249 static const char rmsd_q_doc[] =
01250 "Calculate the root-mean-square distance (RMSD) between selection after\n"
01251 "rotating them optimally\n\n"
01252 "Args:\n"
01253 "    selection (atomsel): Other selection to compute RMSD to. Must have\n"
01254 "        the same number of atoms as this selection\n"
01255 "    weight (list of float): Weight per atom, optional\n"
01256 "Returns:\n"
01257 "    (float): RMSD between selections.";
01258 static PyObject *py_rmsd_q(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01259 {
01260   AtomSel *sel2;
01261   float *weight;
01262   float rmsd;
01263   int rc;
01264 
01265   weight = parse_two_selections_return_weight(a, args, kwargs, &sel2);
01266   if (!weight)
01267     return NULL;
01268 
01269   rc = measure_rmsd_qcp(a->app, a->atomSel, sel2, a->atomSel->selected,
01270                         a->atomSel->coordinates(a->app->moleculeList),
01271                         sel2->coordinates(a->app->moleculeList),
01272                         weight, &rmsd);
01273   delete [] weight;
01274 
01275   if (rc < 0) {
01276     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01277     return NULL;
01278   }
01279   return PyFloat_FromDouble(rmsd);
01280 }
01281 
01282 static char *rmsdmat_q_doc = (char *)
01283   "rmsdmat(sel=this, weight=None, first=0, last=-1, step=1) -> pairwise \n"
01284   "  rmsd over a trajectory, returned as a matrix.\n"
01285   "  Weight must be None or list of same size as selection.\n"
01286   "  Uses QCP algorithm internally, so does not depend on pre-alignment.";
01287 static PyObject *py_rmsdmat_q(PyAtomSelObject *a, PyObject *args, PyObject *kwds) {
01288   int first=0, last=-1, step=1;
01289   PyObject *weightobj = NULL;
01290   static char *kwlist[] = {(char *)"weight",(char *)"first",(char *)"last",(char *)"step", NULL };
01291   if (!PyArg_ParseTupleAndKeywords(args, kwds, "|Oiii", kwlist, 
01292         &weightobj, &first, &last, &step)) { return NULL;}
01293 
01294   float *weight = parse_weight(a->atomSel, weightobj);
01295   if (!weight) return NULL;
01296 
01297   if (last < 0)
01298     last = a->app->molecule_numframes(a->atomSel->molid()) - 1;
01299   int framecount = (last - first + 1) / step;
01300   float *rmsdmat = new float[framecount * framecount];
01301 
01302   // compute the rmsd matrix
01303   int rc = measure_rmsdmat_qcp(a->app, a->atomSel, a->app->moleculeList, 
01304                                a->atomSel->selected, weight, 
01305                                first, last, step, rmsdmat);
01306 
01307   delete [] weight;
01308   if (rc < 0) {
01309     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01310     return NULL;
01311   }
01312 
01313   int i,j;
01314   PyObject* returnlist = PyList_New(framecount);
01315   for (i = 0; i < framecount; i++) {
01316     PyObject* listacross = PyList_New(framecount);
01317     for (j = 0; j < framecount; j++) {
01318       PyList_SetItem(listacross, j, Py_BuildValue("f", rmsdmat[i*framecount + j]));
01319     }
01320     PyList_SET_ITEM(returnlist, i, listacross);
01321   }
01322   delete [] rmsdmat;
01323   return returnlist;
01324 }
01325 
01326 static const char rmsf_doc[] =
01327 "Measures the root-mean-square fluctuation (RMSF) along a trajectory on a\n"
01328 "per-atom basis. RMSF is the mean deviation from the average position\n\n"
01329 "Args:\n"
01330 "    first (int): First frame to include. Defaults to 0 (beginning).\n"
01331 "    last (int): Last frame to include. Defaults to -1 (end).\n"
01332 "    step (int): Use every step-th frame. Defaults to 1 (all frames)\n"
01333 "Returns:\n"
01334 "    (list of float): RMSF for each atom in selection";
01335 static PyObject *py_rmsf(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01336 {
01337   const char *kwlist[] = {"first", "last", "step", NULL};
01338   int first = 0, last = -1, step = 1;
01339   PyObject *returnlist = NULL;
01340   int i, ret_val;
01341   float *rmsf;
01342 
01343   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|iii:atomsel.rmsf",
01344                                    (char**) kwlist, &first, &last, &step))
01345     return NULL;
01346 
01347   // Check molecule is valid
01348   if (!get_molecule(a)) {
01349     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01350     return NULL;
01351   }
01352 
01353   rmsf = new float[a->atomSel->selected];
01354   ret_val = measure_rmsf(a->atomSel, a->app->moleculeList, first, last,
01355                          step, rmsf);
01356   if (ret_val < 0) {
01357     PyErr_SetString(PyExc_RuntimeError, measure_error(ret_val));
01358     goto failure;
01359   }
01360 
01361   returnlist = PyList_New(a->atomSel->selected);
01362   for (i = 0; i < a->atomSel->selected; i++) {
01363     PyList_SET_ITEM(returnlist, i, PyFloat_FromDouble(rmsf[i]));
01364     if (PyErr_Occurred()) {
01365       PyErr_SetString(PyExc_RuntimeError, "cannot build rmsd list");
01366       goto failure;
01367     }
01368   }
01369 
01370   delete [] rmsf;
01371   return returnlist;
01372 
01373 failure:
01374   delete [] rmsf;
01375   Py_XDECREF(returnlist);
01376   return NULL;
01377 }
01378 
01379 static const char rgyr_doc[] =
01380 "Calculate the radius of gyration of this selection\n\n"
01381 "Args:\n"
01382 "    weight (list of float): Per-atom weights to apply during calcuation\n"
01383 "        Must be same size as selection. Optional\n"
01384 "Returns:\n"
01385 "    (float): Radius of gyration";
01386 static PyObject *py_rgyr(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01387 {
01388   const char *kwlist[] = {"weight", NULL};
01389   PyObject *weightobj = NULL;
01390   float *weight;
01391   float rgyr;
01392   int rc;
01393 
01394   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|O:atomsel.rgyr",
01395                                    (char**) kwlist, &weightobj))
01396     return NULL;
01397 
01398   weight = parse_weight(a->atomSel, weightobj);
01399   if (!weight)
01400     return NULL;
01401 
01402   rc = measure_rgyr(a->atomSel, a->app->moleculeList, weight, &rgyr);
01403   delete [] weight;
01404 
01405   if (rc < 0) {
01406     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01407     return NULL;
01408   }
01409 
01410   return PyFloat_FromDouble(rgyr);
01411 }
01412 
01413 static const char fit_doc[] =
01414 "Compute the transformation matrix for the root-mean-square (RMS) alignment\n"
01415 "of this selection to the one given. The format of the matrix is suitable\n"
01416 "for passing to the `atomsel.move()` method\n\n"
01417 "Args:\n"
01418 "    selection (atomsel): Selection to compute fit to. Must have the same\n"
01419 "        number of atoms as this selection\n"
01420 "    weight (list of float): Per-atom weights to apply during calculation\n"
01421 "        Must be the same size as this selection. Optional\n"
01422 "Returns:\n"
01423 "    (16-tuple of float): Transformation matrix, in column major / fortran\n"
01424 "        ordering";
01425 static PyObject *py_fit(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01426 {
01427   PyObject *result;
01428   AtomSel *sel2;
01429   float *weight;
01430   Matrix4 mat;
01431   int rc, i;
01432 
01433   weight = parse_two_selections_return_weight(a, args, kwargs, &sel2);
01434   if (!weight)
01435     return NULL;
01436 
01437   rc = measure_fit(a->atomSel, sel2,
01438                    a->atomSel->coordinates(a->app->moleculeList),
01439                    sel2->coordinates(a->app->moleculeList),
01440                    weight, NULL, &mat);
01441   delete [] weight;
01442 
01443   if (rc < 0) {
01444     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01445     return NULL;
01446   }
01447 
01448   result = PyTuple_New(16);
01449   for (i=0; i<16; i++) {
01450     PyTuple_SET_ITEM(result, i, PyFloat_FromDouble(mat.mat[i]));
01451 
01452     if (PyErr_Occurred()) {
01453       PyErr_SetString(PyExc_RuntimeError, "problem building fit matrix");
01454       Py_XDECREF(result);
01455       return NULL;
01456     }
01457   }
01458 
01459   return result;
01460 }
01461 
01462 static const char moveby_doc[] =
01463 "Shift the selection by a vector\n\n"
01464 "Args:\n"
01465 "    vector (3-tuple of float): (x, y, z) movement to apply";
01466 static PyObject *py_moveby(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01467 {
01468   const char *kwlist[] = {"vector", NULL};
01469   AtomSel *atomSel = a->atomSel;
01470   DrawMolecule *mol;
01471   float offset[3];
01472   float *pos;
01473 
01474   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "(fff):atomsel.moveby",
01475                                    (char**) kwlist, &offset[0], &offset[1],
01476                                    &offset[2]))
01477     return NULL;
01478 
01479   if (!(mol = get_molecule(a))) {
01480     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01481     return NULL;
01482   }
01483 
01484   pos = atomSel->coordinates(a->app->moleculeList);
01485   if (!atomSel->selected || !pos) {
01486     PyErr_Format(PyExc_ValueError, "No coordinates in selection '%s'"
01487                  " molid %d", atomSel->cmdStr, atomSel->molid());
01488     return NULL;
01489   }
01490 
01491   for (int i = atomSel->firstsel; i <= atomSel->lastsel; i++) {
01492     if (atomSel->on[i]) {
01493       float *npos = pos + i * 3;
01494       vec_add(npos, npos, offset);
01495     }
01496   }
01497 
01498   mol->force_recalc(DrawMolItem::MOL_REGEN);
01499   Py_INCREF(Py_None);
01500   return Py_None;
01501 }
01502 
01503 static const char move_doc[] =
01504 "Apply a coordinate transformation to the selection. To undo the move,\n"
01505 "calculate the inverse coordinate transformation matrix with\n"
01506 "`numpy.linalg.inv(matrix)` and pass that to this method\n\n"
01507 "Args:\n"
01508 "    matrix (numpy 16, matrix): Coordinate transformation, in form returned\n"
01509 "        by `atomsel.fit()`, column major / fortran ordering";
01510 static PyObject *py_move(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01511 {
01512   const char *kwlist[] = {"matrix", NULL};
01513   AtomSel *atomSel = a->atomSel;
01514   DrawMolecule *mol;
01515   PyObject *matobj;
01516   Matrix4 mat;
01517   int err;
01518 
01519   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O:atomsel.move",
01520                                    (char**) kwlist, &matobj))
01521     return NULL;
01522 
01523   if (!(mol = get_molecule(a))) {
01524     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01525     return NULL;
01526   }
01527 
01528   if (!atomSel->selected || !atomSel->coordinates(a->app->moleculeList)) {
01529     PyErr_Format(PyExc_ValueError, "No coordinates in selection '%s'"
01530                  " molid %d", atomSel->cmdStr, atomSel->molid());
01531     return NULL;
01532   }
01533 
01534   // Exception set inside function call if an error is returned
01535   if (!py_get_vector(matobj, 16, mat.mat))
01536     return NULL;
01537 
01538   err = measure_move(atomSel, atomSel->coordinates(a->app->moleculeList), mat);
01539   if (err) {
01540     PyErr_SetString(PyExc_ValueError, measure_error(err));
01541     return NULL;
01542   }
01543 
01544   mol->force_recalc(DrawMolItem::MOL_REGEN);
01545   Py_INCREF(Py_None);
01546   return Py_None;
01547 }
01548 
01549 static const char contacts_doc[] =
01550 "Finds all atoms in selection within a given distance of any atom in the\n"
01551 "given selection that are not directly bonded to it. Selections can be in\n"
01552 "different molecules.\n\n"
01553 "Args:\n"
01554 "    selection (atomsel): Atom selection to compare against\n"
01555 "    cutoff (float): Distance cutoff for atoms to be considered contacting\n"
01556 "Returns:\n"
01557 "    (2 lists): Atom indices in this selection, and in given selection\n"
01558 "        that are within the cutoff.";
01559 static PyObject *contacts(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01560 {
01561   const char *kwlist[] = {"selection", "cutoff", NULL};
01562   PyObject *result = NULL, *list1 = NULL, *list2 = NULL, *obj2 = NULL;
01563   GridSearchPair *pairlist, *tmp;
01564   const float *ts1, *ts2;
01565   AtomSel *sel1, *sel2;
01566   DrawMolecule *mol;
01567   float cutoff;
01568 
01569   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!f:atomsel.contacts",
01570                                    (char**) kwlist, &Atomsel_Type, &obj2,
01571                                    &cutoff))
01572     return NULL;
01573 
01574   if (!(mol = get_molecule(a))) {
01575     PyErr_SetString(PyExc_ValueError,
01576                     "this selection is in a deleted molecule");
01577     return NULL;
01578   }
01579 
01580   if (!(get_molecule((PyAtomSelObject *)obj2))) {
01581     PyErr_SetString(PyExc_ValueError,
01582                     "other selection is in a deleted molecule");
01583     return NULL;
01584   }
01585   sel1 = a->atomSel;
01586   sel2 = ((PyAtomSelObject *)obj2)->atomSel;
01587 
01588   if (!(sel1->selected) ||
01589       !(ts1 = sel1->coordinates(a->app->moleculeList))) {
01590     PyErr_Format(PyExc_ValueError, "No coordinates in selection '%s'"
01591                  " molid %d", sel1->cmdStr, sel1->molid());
01592     return NULL;
01593   }
01594 
01595   if (!(sel2->selected) ||
01596       !(ts2 = sel2->coordinates(a->app->moleculeList))) {
01597     PyErr_Format(PyExc_ValueError, "No coordinates in selection '%s'"
01598                  " molid %d", sel2->cmdStr, sel2->molid());
01599     return NULL;
01600   }
01601 
01602   // Check cutoffs are valid
01603   if (cutoff <= 0) {
01604     PyErr_SetString(PyExc_ValueError, "cutoff must be > 0");
01605     return NULL;
01606   }
01607 
01608   pairlist = vmd_gridsearch3(ts1, sel1->num_atoms, sel1->on, ts2,
01609                              sel2->num_atoms, sel2->on, cutoff, -1,
01610                              (sel1->num_atoms + sel2->num_atoms) * 27);
01611 
01612   list1 = PyList_New(0);
01613   list2 = PyList_New(0);
01614   if (PyErr_Occurred())
01615     goto failure;
01616 
01617   for (; pairlist; pairlist = tmp) {
01618     // throw out pairs that are already bonded
01619     MolAtom *a1 = mol->atom(pairlist->ind1);
01620     if (sel1->molid() != sel2->molid() || !a1->bonded(pairlist->ind2)) {
01621       // Since PyList_Append does *not* steal a reference, we need to
01622       // keep the object then decref it after append, to not leak references
01623       PyObject *tmp1 = as_pyint(pairlist->ind1);
01624       PyObject *tmp2 = as_pyint(pairlist->ind2);
01625 
01626       PyList_Append(list1, tmp1);
01627       PyList_Append(list2, tmp2);
01628 
01629       Py_DECREF(tmp1);
01630       Py_DECREF(tmp2);
01631 
01632       if (PyErr_Occurred())
01633         goto failure;
01634     }
01635     tmp = pairlist->next;
01636     free(pairlist);
01637   }
01638 
01639   result = PyList_New(2);
01640   PyList_SET_ITEM(result, 0, list1);
01641   PyList_SET_ITEM(result, 1, list2);
01642   if (PyErr_Occurred())
01643     goto failure;
01644 
01645   return result;
01646 
01647 failure:
01648   PyErr_SetString(PyExc_RuntimeError, "Problem building contacts lists");
01649   Py_XDECREF(list1);
01650   Py_XDECREF(list2);
01651   Py_XDECREF(result);
01652 
01653   // Free linked list of GridSearchPairs, if iteration was broken out of
01654   for (; pairlist; pairlist = tmp) {
01655     tmp = pairlist->next;
01656     free(pairlist);
01657   }
01658   return NULL;
01659 }
01660 
01661 static const char py_hbonds_doc[] =
01662 "Get hydrogen bonds present in current frame of selection using simple\n"
01663 "geometric criteria.\n\n"
01664 "Args:\n"
01665 "    cutoff (float): Distance cutoff between donor and acceptor atoms\n"
01666 "    maxangle (float): Angle cutoff between donor, hydrogen, and acceptor.\n"
01667 "        Angle must be less than this value from 180 degrees.\n"
01668 "    acceptor (atomsel): If given, atomselection for selector atoms, and this\n"
01669 "        selection is assumed to have donor atoms. Both selections must be in\n"
01670 "        the same molecule. If there is overlap between donor and acceptor\n"
01671 "        selection, the output may be inaccurate. Optional.\n"
01672 "Returns:\n"
01673 "    (list of 3 lists): Donor atom indices, acceptor atom indices, and\n"
01674 "        proton atom indices of identified hydrogen bonds\n";
01675 // TODO: make it return a dict
01676 static PyObject *py_hbonds(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01677 {
01678   const char *kwlist[] = {"cutoff", "maxangle", "acceptor", NULL};
01679   PyObject *newdonlist, *newacclist, *newhydlist;
01680   PyObject *obj2 = NULL, *result = NULL;
01681   AtomSel *sel1 = a->atomSel, *sel2;
01682   int *donlist, *hydlist, *acclist;
01683   double cutoff, maxangle;
01684   int maxsize, rc, i;
01685   DrawMolecule *mol;
01686   const float *pos;
01687 
01688   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "dd|O!:atomsel.hbonds",
01689                                    (char**) kwlist, &cutoff, &maxangle,
01690                                    &Atomsel_Type, &obj2))
01691     return NULL;
01692 
01693   if (!(mol = get_molecule(a))) {
01694     PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
01695     return NULL;
01696   }
01697 
01698   // Acceptor must be an atomsel
01699   if (obj2 && !atomsel_Check(obj2)) {
01700     PyErr_SetString(PyExc_TypeError, "acceptor must be an atom selection");
01701     return NULL;
01702   }
01703   sel2 = obj2 ? ((PyAtomSelObject *)obj2)->atomSel : NULL;
01704 
01705   if (sel2 && !get_molecule((PyAtomSelObject*) obj2)) {
01706     PyErr_SetString(PyExc_ValueError, "acceptor selection is on a deleted molecule");
01707     return NULL;
01708   }
01709 
01710   // Selections should be on same molecule
01711   if (obj2 && mol != get_molecule((PyAtomSelObject*) obj2)) {
01712     PyErr_SetString(PyExc_ValueError, "acceptor selection must be in the same "
01713                     "molecule as this selection");
01714     return NULL;
01715   }
01716 
01717   if (!(a->atomSel->selected)
01718    || !(pos = sel1->coordinates(a->app->moleculeList))) {
01719     PyErr_Format(PyExc_ValueError, "No coordinates in selection '%s'"
01720                  " molid %d", a->atomSel->cmdStr, a->atomSel->molid());
01721     return NULL;
01722   }
01723 
01724   // Check cutoffs are valid
01725   if (cutoff <= 0) {
01726     PyErr_SetString(PyExc_ValueError, "cutoff must be > 0");
01727     return NULL;
01728   }
01729 
01730   if (maxangle < 0) {
01731     PyErr_SetString(PyExc_ValueError, "maxangle must be non-negative");
01732     return NULL;
01733   }
01734 
01735   // This heuristic is based on ice, where there are < 2 hydrogen bonds per
01736   // atom if hydrogens are in the selection, and exactly 2 if hydrogens are
01737   // not considered.
01738   maxsize = 2 * sel1->num_atoms;
01739   donlist = new int[maxsize];
01740   hydlist = new int[maxsize];
01741   acclist = new int[maxsize];
01742   rc = measure_hbonds((Molecule *)mol, sel1, sel2, cutoff, maxangle, donlist,
01743                       hydlist, acclist, maxsize);
01744   if (rc > maxsize) {
01745     delete [] donlist;
01746     delete [] hydlist;
01747     delete [] acclist;
01748     maxsize = rc;
01749     donlist = new int[maxsize];
01750     hydlist = new int[maxsize];
01751     acclist = new int[maxsize];
01752     rc = measure_hbonds((Molecule *)mol, sel1, sel2, cutoff, maxangle, donlist,
01753                         hydlist, acclist, maxsize);
01754   }
01755   if (rc < 0) {
01756     PyErr_SetString(PyExc_RuntimeError, "Problem calculating hbonds");
01757     return NULL;
01758   }
01759 
01760   newdonlist = PyList_New(rc);
01761   newacclist = PyList_New(rc);
01762   newhydlist = PyList_New(rc);
01763   if (PyErr_Occurred())
01764     goto failure;
01765 
01766   for (i = 0; i < rc; i++) {
01767     PyList_SET_ITEM(newdonlist, i, as_pyint(donlist[i]));
01768     PyList_SET_ITEM(newacclist, i, as_pyint(acclist[i]));
01769     PyList_SET_ITEM(newhydlist, i, as_pyint(hydlist[i]));
01770 
01771     if (PyErr_Occurred())
01772       goto failure;
01773   }
01774 
01775   result = PyList_New(3);
01776   PyList_SET_ITEM(result, 0, newdonlist);
01777   PyList_SET_ITEM(result, 1, newacclist);
01778   PyList_SET_ITEM(result, 2, newhydlist);
01779 
01780   if (PyErr_Occurred())
01781     goto failure;
01782 
01783   delete [] donlist;
01784   delete [] hydlist;
01785   delete [] acclist;
01786   return result;
01787 
01788 failure:
01789   PyErr_SetString(PyExc_RuntimeError, "Problem building hbonds result");
01790   delete [] donlist;
01791   delete [] hydlist;
01792   delete [] acclist;
01793   Py_XDECREF(newdonlist);
01794   Py_XDECREF(newacclist);
01795   Py_XDECREF(newhydlist);
01796   Py_XDECREF(result);
01797   return NULL;
01798 }
01799 
01800 static const char sasa_doc[] =
01801 "Get solvent accessible surface area of selection\n\n"
01802 "Args:\n"
01803 "    srad (float): Solvent radius\n"
01804 "    samples (int): Maximum number of sample points per atom. Defaults to 500\n"
01805 "    points (bool): True to also return the coordinates of surface points.\n"
01806 "        Defaults to True.\n"
01807 "    restrict (atomsel): Calculate SASA contributions from atoms in this\n"
01808 "        selection only. Useful for getting SASA of residues in the context\n"
01809 "        of their environment. Optional\n"
01810 "Returns:\n"
01811 "    (float): Solvent accessible surface area\n"
01812 "    OR (float, list of 3-tuple): SASA, points, if points=True\n"
01813 "Example to get percent solvent accssibility of a ligand\n"
01814 "    >>> big_sel = atomsel('protein or resname LIG')\n"
01815 "    >>> lig_sel = atomsel('resname LIG')\n"
01816 "    >>> ligand_in_protein_sasa = big_sel.sasa(srad=1.4, restrict=lig_sel)\n"
01817 "    >>> ligand_alone_sasa, points = lig_sel.sasa(srad=1.4, points=True)\n"
01818 "    >>> print(100. * ligand_in_protein_sasa / ligand_alone_sasa)";
01819 static PyObject *sasa(PyAtomSelObject *a, PyObject *args, PyObject *kwargs)
01820 {
01821   const char *kwlist[] = {"srad", "samples", "points", "restrict", NULL};
01822   PyObject *restrictobj = NULL, *pointsobj = NULL;
01823   AtomSel *sel, *restrictsel;
01824   const float *radii, *coords;
01825   ResizeArray<float> sasapts;
01826   float srad = 0, sasa = 0;
01827   int samples = 500, points = 0;
01828   DrawMolecule *mol;
01829   int rc;
01830 
01831   if (!(mol = get_molecule(a)))
01832     return NULL;
01833 
01834   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "f|iO&O!:atomsel.sasa",
01835                                    (char**) kwlist, &srad, &samples,
01836                                    convert_bool, &points, &Atomsel_Type,
01837                                    &restrictobj))
01838     return NULL;
01839 
01840   if (srad < 0) {
01841     PyErr_SetString(PyExc_ValueError, "srad must be non-negative.");
01842     return NULL;
01843   }
01844 
01845   if (samples <= 0) {
01846     PyErr_SetString(PyExc_ValueError, "samples must be positive.");
01847     return NULL;
01848   }
01849 
01850   // fetch the radii and coordinates
01851   sel = a->atomSel;
01852   radii = mol->extraflt.data("radius");
01853   coords = sel->coordinates(a->app->moleculeList);
01854 
01855   // if restrict is given, validate it and pull the selection out
01856   if (restrictobj && !get_molecule((PyAtomSelObject*) restrictobj)) {
01857     PyErr_SetString(PyExc_ValueError, "restrict sel is on deleted molecule");
01858     return NULL;
01859   }
01860   restrictsel = restrictobj ? ((PyAtomSelObject*) restrictobj)->atomSel : NULL;
01861 
01862   // actually calculate sasa
01863   rc = measure_sasa(sel, coords, radii, srad, &sasa,
01864                     points ? &sasapts : NULL, restrictsel, &samples);
01865   if (rc) {
01866     PyErr_SetString(PyExc_ValueError, measure_error(rc));
01867     return NULL;
01868   }
01869 
01870   // append surface points to the provided list object.
01871   if (points) {
01872     pointsobj  = PyList_New(sasapts.num()/3);
01873 
01874     for (int i = 0; i < sasapts.num() / 3; i++) {
01875       PyObject *coord = Py_BuildValue("ddd", sasapts[3L*i], sasapts[3L*i+1],
01876                                       sasapts[3L*i+2]);
01877       PyList_SET_ITEM(pointsobj, i, coord);
01878 
01879       if (PyErr_Occurred()) {
01880         PyErr_SetString(PyExc_RuntimeError, "Problem building points list");
01881         Py_XDECREF(pointsobj);
01882         return NULL;
01883       }
01884     }
01885     return Py_BuildValue("dO", sasa, pointsobj);
01886   }
01887 
01888   return PyFloat_FromDouble(sasa);
01889 }
01890 
01891 #if defined(VMDCUDA)
01892 #include "CUDAMDFF.h"
01893 static const char mdffsim_doc[] =
01894 "Compute a density map\n\n"
01895 "Args:\n"
01896 "    resolution (float): Resolution. Defaults to 10.0\n"
01897 "    spacing (float): Space between adjacent voxel points.\n"
01898 "        Defaults to 0.3*resolution\n"
01899 "Returns:\n"
01900 "    (4-element list): Density map, consisting of 4 lists:\n"
01901 "        1) A 1D list of the grid values at each point\n"
01902 "        2) 3 integers describing the x, y, z lengths in number of grid cells\n"
01903 "        3) 3 floats describing the (x,y,z) coordinates of the origin\n"
01904 "        4) 9 floats describing the deltas for the x, y, and z axes\n"
01905 "Example usage for export to numpy:\n"
01906 "    >>> data, shape, origin, delta = asel.mdffsim(10,3)\n"
01907 "    >>> data = np.array(data)\n"
01908 "    >>> shape = np.array(shape)\n"
01909 "    >>> data = data.reshape(shape, order='F')\n"
01910 "    >>> delta = np.array(delta).reshape(3,3)\n"
01911 "    >>> delta /= shape-1";
01912 static PyObject *py_mdffsim(PyAtomSelObject *a, PyObject *args,
01913                             PyObject *kwargs)
01914 {
01915   const char *kwlist[] = {"resolution", "spacing", NULL};
01916   PyObject *data = NULL, *deltas = NULL, *origin = NULL, *size = NULL;
01917   double gspacing = 0, resolution = 10.0;
01918   VolumetricData *synthvol = NULL;
01919   AtomSel *sel = a->atomSel;
01920   int i, cuda_err, quality;
01921   Molecule *mymol;
01922   float radscale;
01923 
01924   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|dd:atomsel:mdffsim",
01925                                    (char**) kwlist, &resolution, &gspacing))
01926     return NULL;
01927 
01928   if (gspacing < 0) {
01929     PyErr_SetString(PyExc_ValueError, "Spacing must be non-negative");
01930     return NULL;
01931   }
01932 
01933   if (resolution < 0) {
01934     PyErr_SetString(PyExc_ValueError, "Resolution must be non-negative");
01935     return NULL;
01936   }
01937 
01938   if (!(mymol = a->app->moleculeList->mol_from_id(sel->molid()))) {
01939     PyErr_SetString(PyExc_ValueError, "Selection points to a deleted molecule.");
01940     return NULL;
01941   }
01942 
01943   quality = resolution >= 9 ? 0 : 3;
01944   radscale = .2*resolution;
01945   if (gspacing == 0) {
01946     gspacing = 1.5*radscale;
01947   }
01948 
01949   cuda_err = vmd_cuda_calc_density(sel, a->app->moleculeList, quality, radscale,
01950                                    gspacing, &synthvol, NULL, NULL, 1);
01951   if (cuda_err == -1) {
01952     PyErr_SetString(PyExc_ValueError, "CUDA Error, no map returned.");
01953     return NULL;
01954   }
01955 
01956   // after we made it through all of the host code for common error
01957   // handling, we allocate the list object at the point where we can 
01958   // ensure it is properly cleaned up in the error handling branch
01959   PyObject *result = PyList_New(4);
01960 
01961   data = PyList_New(synthvol->xsize * synthvol->ysize * synthvol->zsize);
01962   for (i = 0; i < synthvol->xsize * synthvol->ysize * synthvol->zsize; i++) {
01963     PyList_SET_ITEM(data, i, PyFloat_FromDouble((double)synthvol->data[i]));
01964     if (PyErr_Occurred())
01965       goto failure;
01966   }
01967 
01968   deltas = PyList_New(9);
01969   origin = PyList_New(3);
01970   size = PyList_New(3);
01971   if (PyErr_Occurred())
01972     goto failure;
01973 
01974   // (x, y, z) size
01975   PyList_SET_ITEM(size, 0, as_pyint(synthvol->xsize));
01976   PyList_SET_ITEM(size, 1, as_pyint(synthvol->ysize));
01977   PyList_SET_ITEM(size, 2, as_pyint(synthvol->zsize));
01978   if (PyErr_Occurred())
01979     goto failure;
01980 
01981   // Build length 9 array of axis deltas
01982   for (i = 0; i < 3; i++) {
01983     PyList_SET_ITEM(deltas, i, PyFloat_FromDouble(synthvol->xaxis[i]));
01984     PyList_SET_ITEM(deltas, 3+i, PyFloat_FromDouble(synthvol->yaxis[i]));
01985     PyList_SET_ITEM(deltas, 6+i, PyFloat_FromDouble(synthvol->zaxis[i]));
01986     PyList_SET_ITEM(origin, i, PyFloat_FromDouble(synthvol->origin[i]));
01987     if (PyErr_Occurred())
01988       goto failure;
01989   }
01990 
01991   delete synthvol;
01992   synthvol = NULL;
01993 
01994   PyList_SET_ITEM(result, 0, data);
01995   PyList_SET_ITEM(result, 1, size);
01996   PyList_SET_ITEM(result, 2, origin);
01997   PyList_SET_ITEM(result, 3, deltas);
01998   if (PyErr_Occurred())
01999     goto failure;
02000 
02001   return result;
02002 
02003 failure:
02004   PyErr_SetString(PyExc_RuntimeError, "Problem building grid");
02005   Py_XDECREF(data);
02006   Py_XDECREF(deltas);
02007   Py_XDECREF(origin);
02008   Py_XDECREF(size);
02009   Py_XDECREF(result);
02010 
02011   if(synthvol)
02012     delete synthvol;
02013 
02014   return NULL;
02015 }
02016 
02017 static const char mdffcc_doc[] =
02018 "Get the cross correlation between a volumetric map and the current "
02019 "selection.\n\n"
02020 "Args:\n"
02021 "    volid (int): Index of volumetric dataset\n"
02022 "    resolution (float): Resolution. Defaults to 10.0\n"
02023 "    spacing (float): Space between adjacent voxel points.\n"
02024 "        Defaults to 0.3*resolution\n"
02025 "Returns:\n"
02026 "    (float): Cross correlation for a single frame";
02027 static PyObject *py_mdffcc(PyAtomSelObject *a, PyObject *args,
02028                            PyObject *kwargs)
02029 {
02030   const char *kwlist[] = {"volid", "resolution", "spacing", NULL};
02031   const VolumetricData *volmapA = NULL;
02032   int volid, quality, cuda_err;
02033   const AtomSel *sel = a->atomSel;
02034   float return_cc, radscale;
02035   double resolution = 10.0;
02036   double gspacing = 0;
02037   Molecule *mymol;
02038 
02039   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "i|dd:atomsel.mdffcc",
02040                                    (char**) kwlist, &volid, &resolution,
02041                                    &gspacing))
02042     return NULL;
02043 
02044   if (gspacing < 0) {
02045     PyErr_SetString(PyExc_ValueError, "Spacing must be non-negative");
02046     return NULL;
02047   }
02048 
02049   if (resolution < 0) {
02050     PyErr_SetString(PyExc_ValueError, "Resolution must be non-negative");
02051     return NULL;
02052   }
02053 
02054   if (!(mymol = a->app->moleculeList->mol_from_id(sel->molid()))) {
02055     PyErr_SetString(PyExc_ValueError, "Selection points to a deleted molecule");
02056     return NULL;
02057   }
02058 
02059   if (!(volmapA = mymol->get_volume_data(volid))) {
02060     PyErr_SetString(PyExc_ValueError, "Invalid volume specified. Make sure it's"
02061                     " loaded into the same molecule as the selection");
02062     return NULL;
02063   }
02064 
02065   quality = resolution >= 9 ? 0 : 3;
02066   radscale = .2*resolution;
02067   if (!gspacing) {
02068     gspacing = 1.5*radscale;
02069   }
02070 
02071 
02072   cuda_err = vmd_cuda_compare_sel_refmap(sel, a->app->moleculeList, quality,
02073                                          radscale, gspacing, volmapA, NULL,
02074                                          NULL, NULL, &return_cc, 0.0f, 0);
02075   if (cuda_err == -1) {
02076     PyErr_SetString(PyExc_ValueError, "CUDA Error, no map returned.");
02077     return NULL;
02078   }
02079 
02080   return PyFloat_FromDouble(return_cc);
02081 }
02082 #endif
02083 
02084 // __len__
02085 static Py_ssize_t
02086 atomselection_length( PyObject *a ) {
02087   return ((PyAtomSelObject *)a)->atomSel->selected;
02088 }
02089 
02090 // for integer argument, return True or False if index in in selection
02091 static PyObject *
02092 atomselection_subscript( PyAtomSelObject * a, PyObject * keyobj ) {
02093 #if PY_MAJOR_VERSION >= 3
02094   long ind = PyLong_AsLong(keyobj);
02095 #else
02096   long ind = PyInt_AsLong(keyobj);
02097 #endif
02098   if (ind < 0 && PyErr_Occurred()) return NULL;
02099   PyObject *result = Py_False;
02100   if (ind >= 0 && ind < a->atomSel->num_atoms &&
02101       a->atomSel->on[ind]) {
02102     result = Py_True;
02103   }
02104   Py_INCREF(result);
02105   return result;
02106 }
02107 
02108 static PyMappingMethods atomsel_mapping = {
02109   atomselection_length,
02110   (binaryfunc)atomselection_subscript,
02111   0
02112 };
02113 
02114 /* Methods on selection instances */
02115 static PyMethodDef atomselection_methods[] = {
02116   { "list_attributes", (PyCFunction)py_list_attrs, METH_VARARGS|METH_KEYWORDS, listattrs_doc },
02117   { "get", (PyCFunction)legacy_atomsel_get, METH_VARARGS|METH_KEYWORDS, get_doc  },
02118   { "set", (PyCFunction)legacy_atomsel_set, METH_VARARGS|METH_KEYWORDS, set_doc },
02119   { "update", (PyCFunction)py_update, METH_NOARGS, update_doc },
02120   { "write", (PyCFunction)py_write, METH_VARARGS|METH_KEYWORDS, write_doc },
02121   { "minmax", (PyCFunction)minmax, METH_VARARGS|METH_KEYWORDS, minmax_doc },
02122   { "center", (PyCFunction)center, METH_VARARGS|METH_KEYWORDS, center_doc },
02123   { "rmsd", (PyCFunction)py_rmsd, METH_VARARGS|METH_KEYWORDS, rmsd_doc },
02124   { "rmsdQCP", (PyCFunction)py_rmsd_q, METH_VARARGS|METH_KEYWORDS, rmsd_q_doc },
02125   { "rmsdmat", (PyCFunction)py_rmsdmat_q, METH_VARARGS|METH_KEYWORDS, rmsdmat_q_doc },
02126   { "rmsf", (PyCFunction)py_rmsf, METH_VARARGS|METH_KEYWORDS, rmsf_doc },
02127   { "centerperresidue", (PyCFunction)centerperresidue, METH_VARARGS|METH_KEYWORDS, centerperres_doc },
02128   { "rmsdperresidue", (PyCFunction)py_rmsdperresidue, METH_VARARGS|METH_KEYWORDS, rmsdperres_doc },
02129   { "rmsfperresidue", (PyCFunction)py_rmsfperresidue, METH_VARARGS|METH_KEYWORDS, rmsfperres_doc },
02130   { "rgyr", (PyCFunction)py_rgyr, METH_VARARGS|METH_KEYWORDS, rgyr_doc },
02131   { "fit", (PyCFunction)py_fit, METH_VARARGS|METH_KEYWORDS, fit_doc },
02132   { "move", (PyCFunction)py_move, METH_VARARGS|METH_KEYWORDS, move_doc },
02133   { "moveby", (PyCFunction)py_moveby, METH_VARARGS|METH_KEYWORDS, moveby_doc },
02134   { "contacts", (PyCFunction)contacts, METH_VARARGS|METH_KEYWORDS, contacts_doc },
02135   { "hbonds", (PyCFunction)py_hbonds, METH_VARARGS|METH_KEYWORDS, py_hbonds_doc },
02136   { "sasa", (PyCFunction)sasa, METH_VARARGS|METH_KEYWORDS, sasa_doc },
02137 #if defined(VMDCUDA)
02138   { "mdffsim", (PyCFunction)py_mdffsim, METH_VARARGS|METH_KEYWORDS, mdffsim_doc },
02139   { "mdffcc", (PyCFunction)py_mdffcc, METH_VARARGS|METH_KEYWORDS, mdffcc_doc },
02140 #endif
02141   { NULL, NULL }
02142 };
02143 
02144 // Atom selection iterator
02145 //
02146 typedef struct {
02147   PyObject_HEAD
02148   int index;
02149   PyAtomSelObject * a;
02150 } atomsel_iterobject;
02151 
02152 PyObject *atomsel_iter(PyObject *);
02153 
02154 PyObject *iter_next(atomsel_iterobject *it) {
02155   for ( ; it->index < it->a->atomSel->num_atoms; ++it->index) {
02156     if (it->a->atomSel->on[it->index])
02157       return as_pyint(it->index++);
02158   }
02159   return NULL;
02160 }
02161 
02162 void iter_dealloc(atomsel_iterobject *it) {
02163   Py_XDECREF(it->a);
02164 }
02165 
02166 // Length
02167 PyObject *iter_len(atomsel_iterobject *it) {
02168   return as_pyint(it->a->atomSel->selected);
02169 }
02170 
02171 PyMethodDef iter_methods[] = {
02172   {"__length_hint__", (PyCFunction)iter_len, METH_NOARGS },
02173   {NULL, NULL}
02174 };
02175 
02176 #if PY_MAJOR_VERSION >= 3
02177   PyTypeObject itertype = {
02178     // XXX Clang++ doesn't like this initializer, see PEP 3123:
02179     //     https://www.python.org/dev/peps/pep-3123/
02180     //     See also:
02181     //       http://python3porting.com/cextensions.html#module-initialization
02182     //       https://docs.python.org/3.7/extending/newtypes.html
02183     // PyObject_HEAD_INIT(&PyType_Type)
02184     PyVarObject_HEAD_INIT(NULL, 0)
02185     "atomsel.iterator",
02186     sizeof(atomsel_iterobject), 0, // basic, item size
02187     (destructor)iter_dealloc, // dealloc
02188     0, //tp_print
02189     0, 0, // tp get and setattr
02190     0, // tp_as_async
02191     0, // tp_repr
02192     0, 0, 0, // as number, sequence, mapping
02193     0, 0, 0, // hash, call, str
02194     PyObject_GenericGetAttr, 0, // getattro, setattro
02195     0, // tp_as_buffer
02196     Py_TPFLAGS_DEFAULT, // flags
02197     0, // docstring
02198     0, 0, 0, // traverse, clear, richcompare
02199     0, // tp_weaklistoffset
02200     PyObject_SelfIter, // tp_iter
02201     (iternextfunc)iter_next, // tp_iternext
02202     iter_methods, // tp_methods
02203     0, 0, 0, // members, getset, base
02204   };
02205 
02206   PyTypeObject Atomsel_Type = {
02207     // XXX Clang++ doesn't like this initializer:
02208     //     https://www.python.org/dev/peps/pep-3123/
02209     //     See also:
02210     //       http://python3porting.com/cextensions.html#module-initialization
02211     //       https://docs.python.org/3.7/extending/newtypes.html
02212     // PyObject_HEAD_INIT(0)
02213     PyVarObject_HEAD_INIT(NULL, 0)
02214     "atomsel",
02215     sizeof(PyAtomSelObject), 0, // basic, item size
02216     (destructor)atomsel_dealloc, //dealloc
02217     0, // tp_print
02218     0, 0, // tp get and set attr
02219     0, // tp_as_async
02220     (reprfunc)atomsel_repr, // tp_repr
02221     0, 0, &atomsel_mapping, // as number, sequence, mapping
02222     0, 0, (reprfunc)atomsel_str, // hash, call, str
02223     (getattrofunc) atomsel_getattro, // getattro
02224     (setattrofunc) atomsel_setattro, // setattro
02225     0, // tp_as_buffer
02226     Py_TPFLAGS_DEFAULT|Py_TPFLAGS_BASETYPE, // flags
02227     atomsel_doc, // docstring
02228     0, 0, 0, // traverse, clear, richcompare
02229     0, // tp_weaklistoffset
02230     atomsel_iter, // tp_iter
02231     0, // tp_iternext
02232     atomselection_methods, // tp_methods,
02233     0, atomsel_getset, 0, // members, getset, base
02234     0, 0, 0, // tp_dict, descr_get, descr_set
02235     0, 0, // dictoffset, init
02236     PyType_GenericAlloc, // tp_alloc
02237     atomsel_new, // tp_new
02238     PyObject_Del, // tp_free
02239   };
02240 
02241 #else
02242 
02243   PyTypeObject itertype = {
02244     PyObject_HEAD_INIT(&PyType_Type)
02245     0,
02246     "atomsel.iterator",
02247     sizeof(atomsel_iterobject),
02248     0, // itemsize
02249     /* methods */
02250     (destructor)iter_dealloc,
02251     0,                           /* tp_print */
02252     0,                           /* tp_getattr */
02253     0,                           /* tp_setattr */
02254     0,                           /* tp_compare */
02255     0,                           /* tp_repr */
02256     0,                           /* tp_as_number */
02257     0,                           /* tp_as_sequence */
02258     0,                           /* tp_as_mapping */
02259     0,                           /* tp_hash */
02260     0,                           /* tp_call */
02261     0,                           /* tp_str */
02262     PyObject_GenericGetAttr,     /* tp_getattro */
02263     0,                           /* tp_setattro */
02264     0,                           /* tp_as_buffer */
02265     Py_TPFLAGS_DEFAULT,          /* tp_flags */
02266     atomsel_doc,                 /* tp_doc */
02267     0,                           /* tp_traverse */
02268     0,                           /* tp_clear */
02269     0,                           /* tp_richcompare */
02270     0,                           /* tp_weaklistoffset */
02271     PyObject_SelfIter,           /* tp_iter */
02272     (iternextfunc)iter_next,     /* tp_iternext */
02273     iter_methods,                /* tp_methods */
02274     0,                           /* tp_members */
02275   };
02276 
02277   PyTypeObject Atomsel_Type = {
02278     PyObject_HEAD_INIT(0)        /* Must fill in type value later */
02279     0,                           /* ob_size */
02280     "atomsel.atomsel",           /* tp_name */
02281     sizeof(PyAtomSelObject),     /* tp_basicsize */
02282     0,                           /* tp_itemsize */
02283     (destructor)atomsel_dealloc, /* tp_dealloc */
02284     0,                           /* tp_print */
02285     0,                           /* tp_getattr - getattro used instead */
02286     0,                           /* tp_setattr  - getattro used instead */
02287     0,                           /* tp_compare */
02288     (reprfunc)atomsel_repr,      /* tp_repr */
02289     0,                           /* tp_as_number */
02290     0,                           /* tp_as_sequence */
02291     &atomsel_mapping,            /* tp_as_mapping */
02292     0,                           /* tp_hash */
02293     0,                           /* tp_call */
02294     (reprfunc)atomsel_str,       /* tp_str */
02295     (getattrofunc) atomsel_getattro, //tp_getattro
02296     (setattrofunc) atomsel_setattro, // setattro
02297     0,                           /* tp_as_buffer */
02298     Py_TPFLAGS_DEFAULT|Py_TPFLAGS_BASETYPE|Py_TPFLAGS_HAVE_CLASS, /* tp_flags */
02299     atomsel_doc,                 /* tp_doc */
02300     0,                           /* tp_traverse */
02301     0,                           /* tp_clear */
02302     0,                           /* tp_richcompare */
02303     0,                           /* tp_weaklistoffset */
02304     atomsel_iter,                /* tp_iter */
02305     0,                           /* tp_iternext */
02306     atomselection_methods,       /* tp_methods */
02307     0,                           /* tp_members */
02308     atomsel_getset,              /* tp_getset */
02309     0,                           /* tp_base */
02310     0,                           /* tp_dict */
02311     0,                           /* tp_descr_get */
02312     0,                           /* tp_descr_set */
02313     0,                           /* tp_dictoffset */
02314     0,                           /* tp_init */
02315     PyType_GenericAlloc,         /* tp_alloc */
02316     atomsel_new,                 /* tp_new */
02317     _PyObject_Del        ,       /* tp_free */
02318   };
02319 
02320 #endif
02321 
02322 // tp_iter for atomsel type
02323 PyObject *atomsel_iter(PyObject *self) {
02324   atomsel_iterobject *iter = PyObject_New(atomsel_iterobject, &itertype);
02325   if (!iter)
02326     return NULL;
02327 
02328   Py_INCREF( iter->a = (PyAtomSelObject *)self );
02329   iter->index = 0;
02330   return (PyObject *)iter;
02331 }
02332 
02333 // Atomsel is a type, not a module. So we just initialize the type
02334 // and return it.
02335 PyObject* initatomsel(void) {
02336 #if PY_MAJOR_VERSION >= 3
02337   ((PyObject*)(&Atomsel_Type))->ob_type = &PyType_Type;
02338 #else
02339   Atomsel_Type.ob_type = &PyType_Type;
02340 #endif
02341 
02342   Py_INCREF((PyObject *)&Atomsel_Type);
02343 
02344   if (PyType_Ready(&Atomsel_Type) < 0) {
02345     PyErr_SetString(PyExc_RuntimeError, "Problem initializing atomsel type");
02346     return NULL;
02347   }
02348 
02349   return (PyObject*) &Atomsel_Type;
02350 }
02351 

Generated on Thu Mar 28 02:44:02 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002