00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
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
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
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
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
00084 static void atomsel_dealloc( PyAtomSelObject *obj ) {
00085 delete obj->atomSel;
00086 ((PyObject *)(obj))->ob_type->tp_free((PyObject *)obj);
00087 }
00088
00089
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
00107 static PyObject* atomsel_str(PyAtomSelObject *obj)
00108 {
00109 return as_pystring(obj->atomSel->cmdStr);
00110 }
00111
00112
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
00176
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
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
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
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);
00233 PyList_SET_ITEM(newlist, j++, val);
00234
00235 if (PyErr_Occurred())
00236 goto failure;
00237 }
00238 }
00239 delete [] tmp;
00240
00241
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
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
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
00340 elem = table->fctns.data(i);
00341 if (elem->is_a != SymbolTableElement::KEYWORD
00342 && elem->is_a != SymbolTableElement::SINGLEWORD)
00343 continue;
00344
00345
00346 if (changeable && !table->is_changeable(i))
00347 continue;
00348
00349
00350 if (!table->fctns.name(i) || !strlen(table->fctns.name(i)))
00351 continue;
00352
00353
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
00384 if (!(result = atomsel_get((PyAtomSelObject*) o, attr)))
00385 return NULL;
00386
00387 return result;
00388 }
00389
00390
00391 static PyObject *atomsel_getattro(PyObject *o, PyObject *attr_name)
00392 {
00393
00394 PyObject *tmp;
00395 if (!(tmp = PyObject_GenericGetAttr(o, attr_name))) {
00396
00397
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
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
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
00438 if (!flgs[i])
00439 continue;
00440
00441
00442 p = (char*) list + i*dtype_size;
00443
00444
00445
00446
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
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
00484 if (!(mol = get_molecule(a))) {
00485 PyErr_SetString(PyExc_ValueError, "selection is on a deleted molecule");
00486 return -1;
00487 }
00488
00489
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
00510 atomsel_ctxt context(table, mol, atomSel->which_frame, NULL);
00511
00512
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
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
00535
00536 }
00537
00538 free(list);
00539 list = NULL;
00540
00541
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
00556 free(list);
00557 return 1;
00558 }
00559
00560
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
00583 if (atomsel_set((PyAtomSelObject*) o, attr, val))
00584 return NULL;
00585
00586 Py_INCREF(Py_None);
00587 return Py_None;
00588 }
00589
00590
00591 static int atomsel_setattro(PyObject *o, PyObject *name, PyObject *value)
00592 {
00593
00594
00595
00596
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
00701
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
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);
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
00863
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
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
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
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
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
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
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
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
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
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
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
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
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
01619 MolAtom *a1 = mol->atom(pairlist->ind1);
01620 if (sel1->molid() != sel2->molid() || !a1->bonded(pairlist->ind2)) {
01621
01622
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
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
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
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
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
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
01736
01737
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
01851 sel = a->atomSel;
01852 radii = mol->extraflt.data("radius");
01853 coords = sel->coordinates(a->app->moleculeList);
01854
01855
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
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
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
01957
01958
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
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
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
02085 static Py_ssize_t
02086 atomselection_length( PyObject *a ) {
02087 return ((PyAtomSelObject *)a)->atomSel->selected;
02088 }
02089
02090
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
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
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
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
02179
02180
02181
02182
02183
02184 PyVarObject_HEAD_INIT(NULL, 0)
02185 "atomsel.iterator",
02186 sizeof(atomsel_iterobject), 0,
02187 (destructor)iter_dealloc,
02188 0,
02189 0, 0,
02190 0,
02191 0,
02192 0, 0, 0,
02193 0, 0, 0,
02194 PyObject_GenericGetAttr, 0,
02195 0,
02196 Py_TPFLAGS_DEFAULT,
02197 0,
02198 0, 0, 0,
02199 0,
02200 PyObject_SelfIter,
02201 (iternextfunc)iter_next,
02202 iter_methods,
02203 0, 0, 0,
02204 };
02205
02206 PyTypeObject Atomsel_Type = {
02207
02208
02209
02210
02211
02212
02213 PyVarObject_HEAD_INIT(NULL, 0)
02214 "atomsel",
02215 sizeof(PyAtomSelObject), 0,
02216 (destructor)atomsel_dealloc,
02217 0,
02218 0, 0,
02219 0,
02220 (reprfunc)atomsel_repr,
02221 0, 0, &atomsel_mapping,
02222 0, 0, (reprfunc)atomsel_str,
02223 (getattrofunc) atomsel_getattro,
02224 (setattrofunc) atomsel_setattro,
02225 0,
02226 Py_TPFLAGS_DEFAULT|Py_TPFLAGS_BASETYPE,
02227 atomsel_doc,
02228 0, 0, 0,
02229 0,
02230 atomsel_iter,
02231 0,
02232 atomselection_methods,
02233 0, atomsel_getset, 0,
02234 0, 0, 0,
02235 0, 0,
02236 PyType_GenericAlloc,
02237 atomsel_new,
02238 PyObject_Del,
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,
02249
02250 (destructor)iter_dealloc,
02251 0,
02252 0,
02253 0,
02254 0,
02255 0,
02256 0,
02257 0,
02258 0,
02259 0,
02260 0,
02261 0,
02262 PyObject_GenericGetAttr,
02263 0,
02264 0,
02265 Py_TPFLAGS_DEFAULT,
02266 atomsel_doc,
02267 0,
02268 0,
02269 0,
02270 0,
02271 PyObject_SelfIter,
02272 (iternextfunc)iter_next,
02273 iter_methods,
02274 0,
02275 };
02276
02277 PyTypeObject Atomsel_Type = {
02278 PyObject_HEAD_INIT(0)
02279 0,
02280 "atomsel.atomsel",
02281 sizeof(PyAtomSelObject),
02282 0,
02283 (destructor)atomsel_dealloc,
02284 0,
02285 0,
02286 0,
02287 0,
02288 (reprfunc)atomsel_repr,
02289 0,
02290 0,
02291 &atomsel_mapping,
02292 0,
02293 0,
02294 (reprfunc)atomsel_str,
02295 (getattrofunc) atomsel_getattro,
02296 (setattrofunc) atomsel_setattro,
02297 0,
02298 Py_TPFLAGS_DEFAULT|Py_TPFLAGS_BASETYPE|Py_TPFLAGS_HAVE_CLASS,
02299 atomsel_doc,
02300 0,
02301 0,
02302 0,
02303 0,
02304 atomsel_iter,
02305 0,
02306 atomselection_methods,
02307 0,
02308 atomsel_getset,
02309 0,
02310 0,
02311 0,
02312 0,
02313 0,
02314 0,
02315 PyType_GenericAlloc,
02316 atomsel_new,
02317 _PyObject_Del ,
02318 };
02319
02320 #endif
02321
02322
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
02334
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