AtomSel.C

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 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: AtomSel.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.181 $ $Date: 2022年01月21日 07:58:29 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 * 
00019 * Parse and maintain the data for selecting atoms.
00020 *
00021 ***************************************************************************/
00022 
00023 #include <math.h>
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #include <ctype.h>
00027 #include <string.h>
00028 
00029 #if defined(ARCH_AIX4)
00030 #include <strings.h>
00031 #endif
00032 
00033 #include "Atom.h"
00034 #include "AtomSel.h"
00035 #include "DrawMolecule.h"
00036 #include "MoleculeList.h"
00037 #include "VMDApp.h"
00038 #include "Inform.h"
00039 #include "SymbolTable.h"
00040 #include "ParseTree.h"
00041 #include "JRegex.h"
00042 #include "VolumetricData.h"
00043 #include "PeriodicTable.h"
00044 #include "DrawRingsUtils.h"
00045 
00046 // 'all'
00047 static int atomsel_all(void * ,int, int *) {
00048 return 1;
00049 }
00050 
00051 // 'none'
00052 static int atomsel_none(void *, int num, int *flgs) {
00053 for (int i=num-1; i>=0; i--) {
00054 flgs[i] = FALSE;
00055 }
00056 return 1;
00057 }
00058 
00059 #define generic_atom_data(fctnname, datatype, field) \
00060 static int fctnname(void *v, int num, datatype *data, int *flgs) { \
00061 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00062 for (int i=0; i<num; i++) { \
00063 if (flgs[i]) { \
00064 data[i] = atom_sel_mol->atom(i)->field; \
00065 } \
00066 } \
00067 return 1; \
00068 }
00069 
00070 
00071 #define generic_set_atom_data(fctnname, datatype, fieldtype, field) \
00072 static int fctnname(void *v, int num, datatype *data, int *flgs) { \
00073 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00074 int i; \
00075 for (i=0; i<num; i++) { \
00076 if (flgs[i]) { \
00077 atom_sel_mol->atom(i)->field = (fieldtype) data[i]; \
00078 } \
00079 } \
00080 return 1; \
00081 }
00082 
00083 
00084 // 'name'
00085 static int atomsel_name(void *v, int num, const char **data, int *flgs) {
00086 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00087 for (int i=0; i<num; i++) {
00088 if (flgs[i])
00089 data[i] = (atom_sel_mol->atomNames).name(
00090 atom_sel_mol->atom(i)->nameindex);
00091 }
00092 return 1;
00093 }
00094 
00095 
00096 // 'type'
00097 static int atomsel_type(void *v, int num, const char **data, int *flgs) {
00098 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00099 for (int i=0; i<num; i++) {
00100 if (flgs[i])
00101 data[i] = (atom_sel_mol->atomTypes).name(
00102 atom_sel_mol->atom(i)->typeindex);
00103 }
00104 return 1;
00105 }
00106 
00107 // XXX probably only use this for internal testing
00108 // 'backbonetype'
00109 static int atomsel_backbonetype(void *v, int num, const char **data, int *flgs) {
00110 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00111 for (int i=0; i<num; i++) {
00112 if (flgs[i]) {
00113 switch (atom_sel_mol->atom(i)->atomType) {
00114 case ATOMNORMAL:
00115 data[i] = "normal";
00116 break;
00117 
00118 case ATOMPROTEINBACK:
00119 data[i] = "proteinback";
00120 break;
00121 
00122 case ATOMNUCLEICBACK:
00123 data[i] = "nucleicback";
00124 break;
00125 
00126 case ATOMHYDROGEN:
00127 data[i] = "hydrogen";
00128 break;
00129 
00130 default:
00131 data[i] = "unknown";
00132 break;
00133 }
00134 }
00135 }
00136 return 1;
00137 }
00138 
00139 
00140 // XXX probably only use this for internal testing
00141 // 'residuetype'
00142 static int atomsel_residuetype(void *v, int num, const char **data, int *flgs) {
00143 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00144 for (int i=0; i<num; i++) {
00145 if (flgs[i]) {
00146 switch (atom_sel_mol->atom(i)->residueType) {
00147 case RESNOTHING:
00148 data[i] = "nothing";
00149 break;
00150 
00151 case RESPROTEIN:
00152 data[i] = "protein";
00153 break;
00154 
00155 case RESNUCLEIC:
00156 data[i] = "nucleic";
00157 break;
00158 
00159 case RESWATERS:
00160 data[i] = "waters";
00161 break;
00162 
00163 default:
00164 data[i] = "unknown";
00165 break;
00166 }
00167 }
00168 }
00169 return 1;
00170 }
00171 
00172 
00173 // 'atomicnumber'
00174 generic_atom_data(atomsel_atomicnumber, int, atomicnumber)
00175 
00176 static int atomsel_set_atomicnumber(void *v, int num, int *data, int *flgs) {
00177 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00178 int i;
00179 for (i=0; i<num; i++) {
00180 if (flgs[i]) {
00181 atom_sel_mol->atom(i)->atomicnumber = (int) data[i];
00182 }
00183 }
00184 // when user sets data fields they are marked as valid fields in BaseMolecule
00185 atom_sel_mol->set_dataset_flag(BaseMolecule::ATOMICNUMBER);
00186 return 1;
00187 }
00188 
00189 
00190 // 'element'
00191 static int atomsel_element(void *v, int num, const char **data, int *flgs) {
00192 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00193 for (int i=0; i<num; i++) {
00194 if (flgs[i])
00195 data[i] = get_pte_label(atom_sel_mol->atom(i)->atomicnumber);
00196 }
00197 return 1;
00198 }
00199 static int atomsel_set_element(void *v, int num, const char **data, int *flgs) {
00200 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00201 for (int i=0; i<num; i++) {
00202 if (flgs[i]) {
00203 int idx = get_pte_idx_from_string((const char *)(data[i]));
00204 atom_sel_mol->atom(i)->atomicnumber = idx;
00205 }
00206 }
00207 // when user sets data fields they are marked as valid fields in BaseMolecule
00208 atom_sel_mol->set_dataset_flag(BaseMolecule::ATOMICNUMBER);
00209 return 1;
00210 }
00211 
00212 
00213 // 'index'
00214 static int atomsel_index(void *v, int num, int *data, int *flgs) { 
00215 for (int i=0; i<num; i++) {
00216 if (flgs[i]) {
00217 data[i] = i; // zero-based
00218 }
00219 }
00220 return 1;
00221 }
00222 
00223 
00224 // 'serial' (one-based atom index)
00225 static int atomsel_serial(void *v, int num, int *data, int *flgs) { 
00226 for (int i=0; i<num; i++) {
00227 if (flgs[i]) {
00228 data[i] = i + 1; // one-based
00229 }
00230 }
00231 return 1;
00232 }
00233 
00234 
00235 // 'fragment'
00236 static int atomsel_fragment(void *v, int num, int *data, int *flgs) {
00237 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00238 for (int i=0; i<num; i++) {
00239 if (flgs[i]) {
00240 int residue = atom_sel_mol->atom(i)->uniq_resid;
00241 data[i] = (atom_sel_mol->residue(residue))->fragment;
00242 }
00243 }
00244 return 1;
00245 }
00246 
00247 
00248 // 'numbonds'
00249 generic_atom_data(atomsel_numbonds, int, bonds)
00250 
00251 
00252 // 'residue'
00253 generic_atom_data(atomsel_residue, int, uniq_resid)
00254 
00255 
00256 // 'resname'
00257 static int atomsel_resname(void *v, int num, const char **data, int *flgs) {
00258 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00259 for (int i=0; i<num; i++) {
00260 if (flgs[i])
00261 data[i] = (atom_sel_mol->resNames).name(
00262 atom_sel_mol->atom(i)->resnameindex);
00263 }
00264 return 1;
00265 }
00266 
00267 
00268 // 'altloc'
00269 static int atomsel_altloc(void *v, int num, const char **data, int *flgs) {
00270 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
00271 for (int i=0; i<num; i++) {
00272 if (flgs[i])
00273 data[i] = (atom_sel_mol->altlocNames).name(
00274 atom_sel_mol->atom(i)->altlocindex);
00275 }
00276 return 1;
00277 }
00278 static int atomsel_set_altloc(void *v, int num, const char **data, int *flgs) {
00279 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
00280 NameList<int> &arr = atom_sel_mol->altlocNames;
00281 for (int i=0; i<num; i++) {
00282 if (flgs[i]) {
00283 int ind = arr.add_name((const char *)(data[i]), arr.num());
00284 atom_sel_mol->atom(i)->altlocindex = ind;
00285 }
00286 }
00287 // when user sets data fields they are marked as valid fields in BaseMolecule
00288 atom_sel_mol->set_dataset_flag(BaseMolecule::ALTLOC);
00289 return 1;
00290 }
00291 
00292 
00293 // 'insertion'
00294 generic_atom_data(atomsel_insertion, const char *, insertionstr)
00295 
00296 
00297 // 'chain'
00298 static int atomsel_chain(void *v, int num, const char **data, int *flgs) {
00299 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
00300 for (int i=0; i<num; i++) {
00301 if (flgs[i])
00302 data[i] = (atom_sel_mol->chainNames).name(
00303 atom_sel_mol->atom(i)->chainindex);
00304 }
00305 return 1;
00306 }
00307 
00308 
00309 // 'segname'
00310 static int atomsel_segname(void *v, int num, const char **data, int *flgs) {
00311 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00312 for (int i=0; i<num; i++) {
00313 if (flgs[i])
00314 data[i] = (atom_sel_mol->segNames).name(
00315 atom_sel_mol->atom(i)->segnameindex);
00316 }
00317 return 1;
00318 }
00319 
00320 
00321 // The next few set functions affect the namelists kept in BaseMolecule
00322 static int atomsel_set_name(void *v, int num, const char **data, int *flgs) {
00323 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00324 NameList<int> &arr = atom_sel_mol->atomNames;
00325 for (int i=0; i<num; i++) {
00326 if (flgs[i]) {
00327 int ind = arr.add_name((const char *)(data[i]), arr.num());
00328 atom_sel_mol->atom(i)->nameindex = ind;
00329 }
00330 }
00331 return 1;
00332 }
00333 
00334 static int atomsel_set_type(void *v, int num, const char **data, int *flgs) {
00335 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00336 NameList<int> &arr = atom_sel_mol->atomTypes;
00337 for (int i=0; i<num; i++) {
00338 if (flgs[i]) {
00339 int ind = arr.add_name((const char *)(data[i]), arr.num());
00340 atom_sel_mol->atom(i)->typeindex = ind;
00341 }
00342 }
00343 return 1;
00344 }
00345 
00346 static int atomsel_set_resname(void *v, int num, const char **data, int *flgs) {
00347 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00348 NameList<int> &arr = atom_sel_mol->resNames;
00349 for (int i=0; i<num; i++) {
00350 if (flgs[i]) {
00351 int ind = arr.add_name((const char *)(data[i]), arr.num());
00352 atom_sel_mol->atom(i)->resnameindex = ind;
00353 }
00354 }
00355 return 1;
00356 }
00357 
00358 static int atomsel_set_chain(void *v, int num, const char **data, int *flgs) {
00359 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00360 NameList<int> &arr = atom_sel_mol->chainNames;
00361 for (int i=0; i<num; i++) {
00362 if (flgs[i]) {
00363 int ind = arr.add_name((const char *)(data[i]), arr.num());
00364 atom_sel_mol->atom(i)->chainindex = ind;
00365 }
00366 }
00367 return 1;
00368 }
00369 
00370 static int atomsel_set_segname(void *v, int num, const char **data, int *flgs) {
00371 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00372 NameList<int> &arr = atom_sel_mol->segNames;
00373 for (int i=0; i<num; i++) {
00374 if (flgs[i]) {
00375 int ind = arr.add_name((const char *)(data[i]), arr.num());
00376 atom_sel_mol->atom(i)->segnameindex = ind;
00377 }
00378 }
00379 return 1;
00380 }
00381 
00382 
00383 // 'radius'
00384 static int atomsel_radius(void *v, int num, double *data, int *flgs) {
00385 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00386 float *radius = atom_sel_mol->radius();
00387 for (int i=0; i<num; i++)
00388 if (flgs[i]) data[i] = radius[i];
00389 return 1;
00390 }
00391 static int atomsel_set_radius(void *v, int num, double *data, int *flgs) {
00392 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00393 float *radius = atom_sel_mol->radius();
00394 for (int i=0; i<num; i++)
00395 if (flgs[i]) radius[i] = (float) data[i];
00396 
00397 // when user sets data fields they are marked as valid fields in BaseMolecule
00398 atom_sel_mol->set_dataset_flag(BaseMolecule::RADIUS);
00399 
00400 // Force min/max atom radii to be updated since we've updated the data
00401 atom_sel_mol->set_radii_changed();
00402 
00403 return 1;
00404 }
00405 
00406 
00407 // 'mass'
00408 static int atomsel_mass(void *v, int num, double *data, int *flgs) {
00409 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00410 float *mass = atom_sel_mol->mass();
00411 for (int i=0; i<num; i++)
00412 if (flgs[i]) data[i] = mass[i];
00413 return 1;
00414 }
00415 static int atomsel_set_mass(void *v, int num, double *data, int *flgs) {
00416 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00417 float *mass = atom_sel_mol->mass();
00418 for (int i=0; i<num; i++)
00419 if (flgs[i]) mass[i] = (float) data[i];
00420 // when user sets data fields they are marked as valid fields in BaseMolecule
00421 atom_sel_mol->set_dataset_flag(BaseMolecule::MASS);
00422 return 1;
00423 }
00424 
00425 
00426 // 'charge'
00427 static int atomsel_charge(void *v, int num, double *data, int *flgs) {
00428 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00429 float *charge = atom_sel_mol->charge();
00430 for (int i=0; i<num; i++)
00431 if (flgs[i]) data[i] = charge[i];
00432 return 1;
00433 }
00434 static int atomsel_set_charge(void *v, int num, double *data, int *flgs) {
00435 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00436 float *charge = atom_sel_mol->charge();
00437 for (int i=0; i<num; i++)
00438 if (flgs[i]) charge[i] = (float) data[i];
00439 // when user sets data fields they are marked as valid fields in BaseMolecule
00440 atom_sel_mol->set_dataset_flag(BaseMolecule::CHARGE);
00441 return 1;
00442 }
00443 
00444 
00445 // 'beta'
00446 static int atomsel_beta(void *v, int num, double *data, int *flgs) {
00447 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00448 float *beta = atom_sel_mol->beta();
00449 for (int i=0; i<num; i++)
00450 if (flgs[i]) data[i] = beta[i];
00451 return 1;
00452 }
00453 static int atomsel_set_beta(void *v, int num, double *data, int *flgs) {
00454 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00455 float *beta = atom_sel_mol->beta();
00456 for (int i=0; i<num; i++)
00457 if (flgs[i]) beta[i] = (float) data[i];
00458 // when user sets data fields they are marked as valid fields in BaseMolecule
00459 atom_sel_mol->set_dataset_flag(BaseMolecule::BFACTOR);
00460 return 1;
00461 }
00462 
00463 
00464 // 'occupancy?'
00465 static int atomsel_occupancy(void *v, int num, double *data, int *flgs) {
00466 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00467 float *occupancy = atom_sel_mol->occupancy();
00468 for (int i=0; i<num; i++)
00469 if (flgs[i]) data[i] = occupancy[i];
00470 return 1;
00471 }
00472 static int atomsel_set_occupancy(void *v, int num, double *data, int *flgs) {
00473 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00474 float *occupancy = atom_sel_mol->occupancy();
00475 for (int i=0; i<num; i++)
00476 if (flgs[i]) occupancy[i] = (float) data[i];
00477 // when user sets data fields they are marked as valid fields in BaseMolecule
00478 atom_sel_mol->set_dataset_flag(BaseMolecule::OCCUPANCY);
00479 return 1;
00480 }
00481 
00482 
00483 // 'flags'
00484 static int atomsel_flags(void *v, int num, int *data, int *flgs, int idx) {
00485 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00486 unsigned char *flags = atom_sel_mol->flags();
00487 
00488 printf("atomflags: %d\n", idx);
00489 for (int i=0; i<num; i++) {
00490 if (flgs[i]) {
00491 data[i] = (flags[i] >> idx) & 0x1U;
00492 // printf(" atom[%d] flags: %d masked: %d\n", i, flags[i], data[i]);
00493 }
00494 }
00495 return 1;
00496 }
00497 static int atomsel_set_flags(void *v, int num, int *data, int *flgs, int idx) {
00498 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00499 unsigned char *flags = atom_sel_mol->flags();
00500 for (int i=0; i<num; i++) {
00501 if (flgs[i]) {
00502 unsigned int val = (data[i] != 0);
00503 flags[i] = (flags[i] & (0xff ^ (0x1U << idx))) | (val << idx);
00504 }
00505 }
00506 
00507 // when user sets data fields they are marked as valid fields in BaseMolecule
00508 atom_sel_mol->set_dataset_flag(BaseMolecule::ATOMFLAGS);
00509 return 1;
00510 }
00511 
00512 //
00513 // provide separate callbacks for the different fields due to limitations of
00514 // the existing internal APIs
00515 //
00516 static int atomsel_flag0(void *v, int num, int *data, int *flgs) {
00517 return atomsel_flags(v, num, data, flgs, 0);
00518 }
00519 static int atomsel_set_flag0(void *v, int num, int *data, int *flgs) {
00520 return atomsel_set_flags(v, num, data, flgs, 0);
00521 }
00522 
00523 static int atomsel_flag1(void *v, int num, int *data, int *flgs) {
00524 return atomsel_flags(v, num, data, flgs, 1);
00525 }
00526 static int atomsel_set_flag1(void *v, int num, int *data, int *flgs) {
00527 return atomsel_set_flags(v, num, data, flgs, 1);
00528 }
00529 
00530 static int atomsel_flag2(void *v, int num, int *data, int *flgs) {
00531 return atomsel_flags(v, num, data, flgs, 2);
00532 }
00533 static int atomsel_set_flag2(void *v, int num, int *data, int *flgs) {
00534 return atomsel_set_flags(v, num, data, flgs, 2);
00535 }
00536 
00537 static int atomsel_flag3(void *v, int num, int *data, int *flgs) {
00538 return atomsel_flags(v, num, data, flgs, 3);
00539 }
00540 static int atomsel_set_flag3(void *v, int num, int *data, int *flgs) {
00541 return atomsel_set_flags(v, num, data, flgs, 3);
00542 }
00543 
00544 static int atomsel_flag4(void *v, int num, int *data, int *flgs) {
00545 return atomsel_flags(v, num, data, flgs, 4);
00546 }
00547 static int atomsel_set_flag4(void *v, int num, int *data, int *flgs) {
00548 return atomsel_set_flags(v, num, data, flgs, 4);
00549 }
00550 
00551 static int atomsel_flag5(void *v, int num, int *data, int *flgs) {
00552 return atomsel_flags(v, num, data, flgs, 5);
00553 }
00554 static int atomsel_set_flag5(void *v, int num, int *data, int *flgs) {
00555 return atomsel_set_flags(v, num, data, flgs, 5);
00556 }
00557 
00558 static int atomsel_flag6(void *v, int num, int *data, int *flgs) {
00559 return atomsel_flags(v, num, data, flgs, 6);
00560 }
00561 static int atomsel_set_flag6(void *v, int num, int *data, int *flgs) {
00562 return atomsel_set_flags(v, num, data, flgs, 6);
00563 }
00564 
00565 static int atomsel_flag7(void *v, int num, int *data, int *flgs) {
00566 return atomsel_flags(v, num, data, flgs, 7);
00567 }
00568 static int atomsel_set_flag7(void *v, int num, int *data, int *flgs) {
00569 return atomsel_set_flags(v, num, data, flgs, 7);
00570 }
00571 
00572 
00573 
00574 // 'resid'
00575 static int atomsel_resid(void *v, int num, int *data, int *flgs) {
00576 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00577 for (int i=0; i<num; i++) {
00578 if (flgs[i]) {
00579 data[i] = atom_sel_mol->atom(i)->resid;
00580 }
00581 }
00582 return 1;
00583 }
00584 static int atomsel_set_resid(void *v, int num, int *data, int *flgs) {
00585 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00586 for (int i=0; i<num; i++) {
00587 if (flgs[i])
00588 atom_sel_mol->atom(i)->resid = data[i];
00589 }
00590 return 1;
00591 }
00592 
00593 
00594 
00595 #define generic_atom_boolean(fctnname, comparison) \
00596 static int fctnname(void *v, int num, int *flgs) { \
00597 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00598 for (int i=0; i<num; i++) { \
00599 if (flgs[i]) { \
00600 flgs[i] = atom_sel_mol->atom(i)->comparison; \
00601 } \
00602 } \
00603 return 1; \
00604 }
00605 
00606 // 'backbone'
00607 static int atomsel_backbone(void *v, int num, int *flgs) {
00608 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00609 for (int i=0; i<num; i++) {
00610 if (flgs[i]) {
00611 const MolAtom *a = atom_sel_mol->atom(i);
00612 flgs[i] = (a->atomType == ATOMPROTEINBACK ||
00613 a->atomType == ATOMNUCLEICBACK);
00614 }
00615 }
00616 return 1;
00617 }
00618 
00619 // 'h' (hydrogen)
00620 generic_atom_boolean(atomsel_hydrogen, atomType == ATOMHYDROGEN)
00621 
00622 
00623 // 'protein'
00624 generic_atom_boolean(atomsel_protein, residueType == RESPROTEIN)
00625 
00626 
00627 // 'nucleic'
00628 generic_atom_boolean(atomsel_nucleic, residueType == RESNUCLEIC)
00629 
00630 
00631 // 'water'
00632 generic_atom_boolean(atomsel_water, residueType == RESWATERS)
00633 
00634 #define generic_sstruct_boolean(fctnname, comparison) \
00635 static int fctnname(void *v, int num, int *flgs) { \
00636 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00637 atom_sel_mol->need_secondary_structure(1); \
00638 for (int i=0; i<num; i++) { \
00639 int s; \
00640 if (flgs[i]) { \
00641 s = atom_sel_mol->residue( \
00642 atom_sel_mol->atom(i)->uniq_resid \
00643 )->sstruct; \
00644 if (!comparison) { \
00645 flgs[i] = 0; \
00646 } \
00647 } \
00648 } \
00649 return 1; \
00650 }
00651 
00652 // once I define a structure, I don't need to recompute; hence the 
00653 // need_secondary_structure(0);
00654 #define generic_set_sstruct_boolean(fctnname, newval) \
00655 static int fctnname(void *v, int num, int *data, int *flgs) { \
00656 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00657 atom_sel_mol->need_secondary_structure(0); \
00658 for (int i=0; i<num; i++) { \
00659 if (flgs[i] && data[i]) { \
00660 atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid) \
00661 ->sstruct = newval; \
00662 } \
00663 } \
00664 return 1; \
00665 }
00666 
00667 
00668 // XXX recursive routines should be replaced by an iterative version with
00669 // it's own stack, so that huge molecules don't overflow the stack
00670 static void recursive_find_sidechain_atoms(BaseMolecule *mol, int *sidechain,
00671 int atom_index) {
00672 // Have I been here before?
00673 if (sidechain[atom_index] == 2) 
00674 return;
00675 
00676 // Is this a backbone atom
00677 MolAtom *atom = mol->atom(atom_index);
00678 if (atom->atomType == ATOMPROTEINBACK ||
00679 atom->atomType == ATOMNUCLEICBACK) return;
00680 
00681 // mark this atom
00682 sidechain[atom_index] = 2;
00683 
00684 // try the atoms to which this is bonded
00685 for (int i=0; i<atom->bonds; i++) {
00686 recursive_find_sidechain_atoms(mol, sidechain, atom->bondTo[i]);
00687 }
00688 }
00689 
00690 // give an array where 1 indicates an atom on the sidechain, find the
00691 // connected atoms which are also on the sidechain
00692 static void find_sidechain_atoms(BaseMolecule *mol, int *sidechain) {
00693 for (int i=0; i<mol->nAtoms; i++) {
00694 if (sidechain[i] == 1) {
00695 recursive_find_sidechain_atoms(mol, sidechain, i);
00696 }
00697 }
00698 }
00699 
00700 
00701 // 'sidechain' is tricky. I start from a protein CA, pick a bond which
00702 // isn't along the backbone (and discard it once if it is a hydrogen),
00703 // and follow the atoms until I stop at another backbone atom or run out
00704 // of atoms
00705 static int atomsel_sidechain(void *v, int num, int *flgs) {
00706 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
00707 const float *mass = atom_sel_mol->mass();
00708 int i;
00709 
00710 // generate a list of the "CB" atoms (or whatever they are)
00711 int *seed = new int[num];
00712 memset(seed, 0, num * sizeof(int));
00713 
00714 // get the CA and HA2 name index
00715 int CA = atom_sel_mol->atomNames.typecode((char *) "CA");
00716 
00717 // for each protein
00718 for (int pfrag=atom_sel_mol->pfragList.num()-1; pfrag>=0; pfrag--) {
00719 // get a residue
00720 Fragment &frag = *(atom_sel_mol->pfragList[pfrag]);
00721 for (int res = frag.num()-1; res >=0; res--) {
00722 // find the CA
00723 int atomid = atom_sel_mol->find_atom_in_residue(CA, frag[res]);
00724 if (atomid < 0) {
00725 msgErr << "atomselection: sidechain: cannot find CA" << sendmsg;
00726 continue;
00727 }
00728 // find at most two neighbors which are not on the backbone
00729 MolAtom *atom = atom_sel_mol->atom(atomid);
00730 int b1 = -1, b2 = -1;
00731 for (i=0; i<atom->bonds; i++) {
00732 int bondtype = atom_sel_mol->atom(atom->bondTo[i])->atomType;
00733 if (bondtype == ATOMNORMAL || bondtype == ATOMHYDROGEN) {
00734 if (b1 == -1) {
00735 b1 = atom->bondTo[i];
00736 } else {
00737 if (b2 == -1) {
00738 b2 = atom->bondTo[i];
00739 } else {
00740 msgErr << "atomselection: sidechain: protein residue index " 
00741 << res << ", C-alpha index " << i << " has more than "
00742 << "two non-backbone bonds; ignoring the others"
00743 << sendmsg;
00744 }
00745 }
00746 }
00747 }
00748 if (b1 == -1) 
00749 continue;
00750 
00751 if (b2 != -1) { // find the right one
00752 // first check the number of atoms and see if I have a lone H
00753 int c1 = atom_sel_mol->atom(b1)->bonds;
00754 int c2 = atom_sel_mol->atom(b2)->bonds;
00755 if (c1 == 1 && c2 > 1) {
00756 b1 = b2;
00757 } else if (c2 == 1 && c1 > 1) {
00758 #if 1
00759 // XXX get rid of bogus self-assignment
00760 seed[b1] = 1; // found the right one; it is b1.
00761 continue; // b1 remains b1
00762 #else
00763 b1 = b1;
00764 #endif
00765 } else if (c1 ==1 && c2 == 1) {
00766 // check the masses
00767 float m1 = mass[b1];
00768 float m2 = mass[b2];
00769 if (m1 > 2.3 && m2 <= 2.3) {
00770 b1 = b2;
00771 } else if (m2 > 2.3 && m1 <= 2.3) {
00772 #if 1
00773 // XXX get rid of bogus self-assignment
00774 seed[b1] = 1; // found the right one; it is b1.
00775 continue; // b1 remains b1
00776 #else
00777 b1 = b1;
00778 #endif
00779 } else if (m1 <= 2.0 && m2 <= 2.3) {
00780 // should have two H's, find the "first" of these
00781 if (strcmp(
00782 (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b1)->nameindex),
00783 (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b2)->nameindex)
00784 ) > 0) {
00785 b1 = b2;
00786 } // else b1 = b1
00787 } else {
00788 msgErr << "atomselect: sidechain: protein residue index " 
00789 << res << ", C-alpha index " << i << " has sidechain-like "
00790 << "atom (indices " << b1 << " and " << b2 << "), and "
00791 << "I cannot determine which to call a sidechain -- "
00792 << "I'll guess" << sendmsg;
00793 if (strcmp(
00794 (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b1)->nameindex),
00795 (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b2)->nameindex)
00796 ) > 0) {
00797 b1 = b2;
00798 }
00799 } // punted
00800 } // checked connections and masses
00801 } // found the right one; it is b1.
00802 seed[b1] = 1;
00803 
00804 } // loop over residues
00805 } // loop over protein fragments
00806 
00807 // do the search for all the sidechain atoms (returned in seed)
00808 find_sidechain_atoms(atom_sel_mol, seed);
00809 
00810 // set the return values
00811 for (i=0; i<num; i++) {
00812 if (flgs[i]) 
00813 flgs[i] = (seed[i] != 0);
00814 }
00815 
00816 delete [] seed;
00817 
00818 return 1;
00819 }
00820 
00821 
00822 // which of these are helices?
00823 generic_sstruct_boolean(atomsel_helix, (s == SS_HELIX_ALPHA ||
00824 s == SS_HELIX_3_10 ||
00825 s == SS_HELIX_PI))
00826 generic_sstruct_boolean(atomsel_alpha_helix, (s == SS_HELIX_ALPHA))
00827 generic_sstruct_boolean(atomsel_3_10_helix, (s == SS_HELIX_3_10))
00828 generic_sstruct_boolean(atomsel_pi_helix, (s == SS_HELIX_PI))
00829 
00830 // Makes the residue into a HELIX_ALPHA
00831 generic_set_sstruct_boolean(atomsel_set_helix, SS_HELIX_ALPHA)
00832 // Makes the residue into a HELIX_3_10
00833 generic_set_sstruct_boolean(atomsel_set_3_10_helix, SS_HELIX_3_10)
00834 // Makes the residue into a HELIX_PI
00835 generic_set_sstruct_boolean(atomsel_set_pi_helix, SS_HELIX_PI)
00836 
00837 // which of these are beta sheets?
00838 generic_sstruct_boolean(atomsel_sheet, (s == SS_BETA ||
00839 s == SS_BRIDGE ))
00840 generic_sstruct_boolean(atomsel_extended_sheet, (s == SS_BETA))
00841 generic_sstruct_boolean(atomsel_bridge_sheet, (s == SS_BRIDGE ))
00842 
00843 // Makes the residue into a BETA
00844 generic_set_sstruct_boolean(atomsel_set_sheet, SS_BETA)
00845 generic_set_sstruct_boolean(atomsel_set_bridge_sheet, SS_BRIDGE)
00846 
00847 // which of these are coils?
00848 generic_sstruct_boolean(atomsel_coil, (s == SS_COIL))
00849 
00850 // Makes the residue into COIL
00851 generic_set_sstruct_boolean(atomsel_set_coil, SS_COIL)
00852 
00853 // which of these are TURNS?
00854 generic_sstruct_boolean(atomsel_turn, (s == SS_TURN))
00855 
00856 // Makes the residue into TURN
00857 generic_set_sstruct_boolean(atomsel_set_turn, SS_TURN)
00858 
00859 // return the sstruct information as a 1 character string
00860 static int atomsel_sstruct(void *v, int num, const char **data, int *flgs) {
00861 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00862 atom_sel_mol->need_secondary_structure(1);
00863 for (int i=0; i<num; i++) {
00864 if (flgs[i]) {
00865 switch(atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid)->sstruct) {
00866 case SS_HELIX_ALPHA: data[i] = "H"; break;
00867 case SS_HELIX_3_10 : data[i] = "G"; break;
00868 case SS_HELIX_PI : data[i] = "I"; break;
00869 case SS_BETA : data[i] = "E"; break;
00870 case SS_BRIDGE : data[i] = "B"; break;
00871 case SS_TURN : data[i] = "T"; break;
00872 default:
00873 case SS_COIL : data[i] = "C"; break;
00874 }
00875 }
00876 }
00877 return 1;
00878 }
00879 
00880 
00881 // set the secondary structure based on a string value
00882 static int atomsel_set_sstruct(void *v, int num, const char **data, int *flgs) {
00883 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
00884 char c;
00885 // Defining it here so remind myself that I don't need STRIDE (or
00886 // whoever) to do it automatically
00887 atom_sel_mol->need_secondary_structure(0);
00888 for (int i=0; i<num; i++) {
00889 if (flgs[i]) {
00890 if (strlen(data[i]) == 0) {
00891 msgErr << "cannot set a 0 length secondary structure string"
00892 << sendmsg;
00893 } else {
00894 c = ((const char *) data[i])[0];
00895 if (strlen(data[i]) > 1) {
00896 while (1) {
00897 if (!strcasecmp((const char *) data[i], "helix")) { c = 'H'; break;}
00898 if (!strcasecmp((const char *) data[i], "alpha")) { c = 'H'; break;}
00899 if (!strcasecmp((const char *) data[i], "alpha_helix")) { c = 'H'; break;}
00900 if (!strcasecmp((const char *) data[i], "alphahelix")) { c = 'H'; break;}
00901 if (!strcasecmp((const char *) data[i], "alpha helix")) { c = 'H'; break;}
00902 
00903 if (!strcasecmp((const char *) data[i], "pi")) { c = 'I'; break;}
00904 if (!strcasecmp((const char *) data[i], "pi_helix")) { c = 'I'; break;}
00905 if (!strcasecmp((const char *) data[i], "pihelix")) { c = 'I'; break;}
00906 if (!strcasecmp((const char *) data[i], "pi helix")) { c = 'I'; break;}
00907 
00908 if (!strcasecmp((const char *) data[i], "310")) { c = 'G'; break;}
00909 if (!strcasecmp((const char *) data[i], "310_helix")) { c = 'G'; break;}
00910 if (!strcasecmp((const char *) data[i], "3_10")) { c = 'G'; break;}
00911 if (!strcasecmp((const char *) data[i], "3")) { c = 'G'; break;}
00912 if (!strcasecmp((const char *) data[i], "310 helix")) { c = 'G'; break;}
00913 if (!strcasecmp((const char *) data[i], "3_10_helix")){ c = 'G'; break;}
00914 if (!strcasecmp((const char *) data[i], "3_10 helix")){ c = 'G'; break;}
00915 if (!strcasecmp((const char *) data[i], "3 10 helix")){ c = 'G'; break;}
00916 if (!strcasecmp((const char *) data[i], "helix_3_10")){ c = 'G'; break;}
00917 if (!strcasecmp((const char *) data[i], "helix 3 10")){ c = 'G'; break;}
00918 if (!strcasecmp((const char *) data[i], "helix3_10")) { c = 'G'; break;}
00919 if (!strcasecmp((const char *) data[i], "helix310")) { c = 'G'; break;}
00920 
00921 if (!strcasecmp((const char *) data[i], "beta")) { c = 'E'; break;}
00922 if (!strcasecmp((const char *) data[i], "betasheet")) { c = 'E'; break;}
00923 if (!strcasecmp((const char *) data[i], "beta_sheet")){ c = 'E'; break;}
00924 if (!strcasecmp((const char *) data[i], "beta sheet")){ c = 'E'; break;}
00925 if (!strcasecmp((const char *) data[i], "sheet")) { c = 'E'; break;}
00926 if (!strcasecmp((const char *) data[i], "strand")) { c = 'E'; break;}
00927 if (!strcasecmp((const char *) data[i], "beta_strand")) { c = 'E'; break;}
00928 if (!strcasecmp((const char *) data[i], "beta strand")) { c = 'E'; break;}
00929 
00930 if (!strcasecmp((const char *) data[i], "turn")) { c = 'T'; break;}
00931 
00932 if (!strcasecmp((const char *) data[i], "coil")) { c = 'C'; break;}
00933 if (!strcasecmp((const char *) data[i], "unknown")) { c = 'C'; break;}
00934 c = 'X';
00935 break;
00936 } // while (1)
00937 }
00938 // and set the value
00939 int s = SS_COIL;
00940 switch ( c ) {
00941 case 'H':
00942 case 'h': s = SS_HELIX_ALPHA; break;
00943 case '3':
00944 case 'G':
00945 case 'g': s = SS_HELIX_3_10; break;
00946 case 'P': // so you can say 'pi'
00947 case 'p':
00948 case 'I':
00949 case 'i': s = SS_HELIX_PI; break;
00950 case 'S': // for "sheet"
00951 case 's':
00952 case 'E':
00953 case 'e': s = SS_BETA; break;
00954 case 'B':
00955 case 'b': s = SS_BRIDGE; break;
00956 case 'T':
00957 case 't': s = SS_TURN; break;
00958 case 'L': // L is from PHD and may be an alternate form
00959 case 'l':
00960 case 'C':
00961 case 'c': s = SS_COIL; break;
00962 default: {
00963 msgErr << "Unknown sstruct assignment: '"
00964 << (const char *) data[i] << "'" << sendmsg;
00965 s = SS_COIL; break;
00966 }
00967 }
00968 atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid)->sstruct = s;
00969 }
00970 }
00971 }
00972 return 1;
00973 }
00974 
00975 
00977 // and leave the rest alone.
00978 // It is slower this way, but easier to understand
00979 static void mark_atoms_given_residue(DrawMolecule *mol, int residue, int *on)
00980 {
00981 ResizeArray<int> *atoms = &(mol->residueList[residue]->atoms);
00982 for (int i= atoms->num()-1; i>=0; i--) {
00983 on[(*atoms)[i]] = TRUE;
00984 }
00985 }
00986 
00987 
00988 // macro for either protein or nucleic fragments
00989 #define fragment_data(fctn, fragtype) \
00990 static int fctn(void *v, int num, int *data, int *) \
00991 { \
00992 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
00993 int *tmp = new int[num]; \
00994 int i; \
00995 for (i=num-1; i>=0; i--) { /* clear the arrays */ \
00996 tmp[i] = 0; \
00997 data[i] = -1; /* default is -1 for 'not a [np]frag' */ \
00998 } \
00999 /* for each fragment */ \
01000 for ( i=atom_sel_mol->fragtype.num()-1; i>=0; i--) { \
01001 /* for each residues of the fragment */ \
01002 int j; \
01003 for (j=atom_sel_mol->fragtype[i]->num()-1; j>=0; j--) { \
01004 /* mark the atoms in the fragment */ \
01005 mark_atoms_given_residue(atom_sel_mol,(*atom_sel_mol->fragtype[i])[j], tmp); \
01006 } \
01007 /* and label them with the appropriate number */ \
01008 for (j=num-1; j>=0; j--) { \
01009 if (tmp[j]) { \
01010 data[j] = i; \
01011 tmp[j] = 0; \
01012 } \
01013 } \
01014 } \
01015 delete [] tmp; \
01016 return 1; \
01017 }
01018 
01019 fragment_data(atomsel_pfrag, pfragList)
01020 fragment_data(atomsel_nfrag, nfragList)
01021 
01022 static Timestep *selframe(DrawMolecule *atom_sel_mol, int which_frame) {
01023 Timestep *r;
01024 switch (which_frame) {
01025 case AtomSel::TS_LAST: r = atom_sel_mol->get_last_frame(); break;
01026 
01027 case AtomSel::TS_NOW : r = atom_sel_mol->current(); break;
01028 default: {
01029 if (!atom_sel_mol->get_frame(which_frame)) {
01030 r = atom_sel_mol->get_last_frame();
01031 
01032 } else {
01033 r = atom_sel_mol->get_frame(which_frame);
01034 }
01035 }
01036 }
01037 return r;
01038 }
01039 
01040 
01041 #if defined(VMDWITHCARBS)
01042 // 'pucker'
01043 static int atomsel_pucker(void *v, int num, double *data, int *flgs) {
01044 int i;
01045 SmallRing *ring;
01046 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
01047 memset(data, 0, num*sizeof(double));
01048 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01049 Timestep *ts = selframe(atom_sel_mol, which_frame);
01050 if (!ts || !ts->pos) {
01051 return 1;
01052 }
01053 
01054 // XXX We're hijacking the ring list in BaseMolecule at present.
01055 // It might be better to build our own independent one, but
01056 // this way there's only one ring list in memory at a time.
01057 atom_sel_mol->find_small_rings_and_links(5, 6);
01058 
01059 for (i=0; i < atom_sel_mol->smallringList.num(); i++) {
01060 ring = atom_sel_mol->smallringList[i];
01061 int N = ring->num();
01062 
01063 // accept rings of size 5 or 6
01064 if (N >= 5 && N <= 6) {
01065 float pucker = hill_reilly_ring_pucker((*ring), ts->pos);
01066 
01067 int j;
01068 for (j=0; j<N; j++) {
01069 int ind = (*ring)[j];
01070 if (flgs[ind])
01071 data[ind] = pucker;
01072 }
01073 }
01074 }
01075 
01076 return 1;
01077 }
01078 #endif
01079 
01080 static int atomsel_user(void *v, int num, double *data, int *flgs) {
01081 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01082 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01083 Timestep *ts = selframe(atom_sel_mol, which_frame);
01084 if (!ts || !ts->user) {
01085 memset(data, 0, num*sizeof(double));
01086 return 1; 
01087 }
01088 for (int i=0; i<num; i++) { 
01089 if (flgs[i]) { 
01090 data[i] = ts->user[i]; 
01091 } 
01092 } 
01093 return 1; 
01094 }
01095 static int atomsel_set_user(void *v, int num, double *data, int *flgs) {
01096 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01097 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01098 Timestep *ts = selframe(atom_sel_mol, which_frame);
01099 if (!ts) return 0;
01100 if (!ts->user) {
01101 ts->user= new float[num];
01102 memset(ts->user, 0, num*sizeof(float));
01103 }
01104 for (int i=0; i<num; i++) { 
01105 if (flgs[i]) { 
01106 ts->user[i] = (float)data[i]; 
01107 } 
01108 } 
01109 return 1; 
01110 }
01111 
01112 static int atomsel_user2(void *v, int num, double *data, int *flgs) {
01113 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01114 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01115 Timestep *ts = selframe(atom_sel_mol, which_frame);
01116 if (!ts || !ts->user2) {
01117 memset(data, 0, num*sizeof(double));
01118 return 1; 
01119 }
01120 for (int i=0; i<num; i++) { 
01121 if (flgs[i]) { 
01122 data[i] = ts->user2[i]; 
01123 } 
01124 } 
01125 return 1; 
01126 }
01127 static int atomsel_set_user2(void *v, int num, double *data, int *flgs) {
01128 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01129 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01130 Timestep *ts = selframe(atom_sel_mol, which_frame);
01131 if (!ts) return 0;
01132 if (!ts->user2) {
01133 ts->user2= new float[num];
01134 memset(ts->user2, 0, num*sizeof(float));
01135 }
01136 for (int i=0; i<num; i++) { 
01137 if (flgs[i]) { 
01138 ts->user2[i] = (float)data[i]; 
01139 } 
01140 } 
01141 return 1; 
01142 }
01143 
01144 static int atomsel_user3(void *v, int num, double *data, int *flgs) {
01145 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01146 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01147 Timestep *ts = selframe(atom_sel_mol, which_frame);
01148 if (!ts || !ts->user3) {
01149 memset(data, 0, num*sizeof(double));
01150 return 1; 
01151 }
01152 for (int i=0; i<num; i++) { 
01153 if (flgs[i]) { 
01154 data[i] = ts->user3[i]; 
01155 } 
01156 } 
01157 return 1; 
01158 }
01159 static int atomsel_set_user3(void *v, int num, double *data, int *flgs) {
01160 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01161 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01162 Timestep *ts = selframe(atom_sel_mol, which_frame);
01163 if (!ts) return 0;
01164 if (!ts->user3) {
01165 ts->user3= new float[num];
01166 memset(ts->user3, 0, num*sizeof(float));
01167 }
01168 for (int i=0; i<num; i++) { 
01169 if (flgs[i]) { 
01170 ts->user3[i] = (float)data[i]; 
01171 } 
01172 } 
01173 return 1; 
01174 }
01175 
01176 static int atomsel_user4(void *v, int num, double *data, int *flgs) {
01177 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01178 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01179 Timestep *ts = selframe(atom_sel_mol, which_frame);
01180 if (!ts || !ts->user4) {
01181 memset(data, 0, num*sizeof(double));
01182 return 1; 
01183 }
01184 for (int i=0; i<num; i++) { 
01185 if (flgs[i]) { 
01186 data[i] = ts->user4[i]; 
01187 } 
01188 } 
01189 return 1; 
01190 }
01191 static int atomsel_set_user4(void *v, int num, double *data, int *flgs) {
01192 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01193 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01194 Timestep *ts = selframe(atom_sel_mol, which_frame);
01195 if (!ts) return 0;
01196 if (!ts->user4) {
01197 ts->user4= new float[num];
01198 memset(ts->user4, 0, num*sizeof(float));
01199 }
01200 for (int i=0; i<num; i++) { 
01201 if (flgs[i]) { 
01202 ts->user4[i] = (float)data[i]; 
01203 } 
01204 } 
01205 return 1; 
01206 }
01207 
01208 
01209 
01210 static int atomsel_xpos(void *v, int num, double *data, int *flgs) {
01211 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01212 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01213 Timestep *ts = selframe(atom_sel_mol, which_frame);
01214 int i;
01215 if (!ts) {
01216 for (i=0; i<num; i++) 
01217 if (flgs[i]) data[i] = 0.0;
01218 return 0; 
01219 }
01220 for (i=0; i<num; i++) { 
01221 if (flgs[i]) { 
01222 data[i] = ts->pos[3L*i ]; 
01223 } 
01224 } 
01225 return 1; 
01226 }
01227 
01228 static int atomsel_ypos(void *v, int num, double *data, int *flgs) {
01229 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01230 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01231 Timestep *ts = selframe(atom_sel_mol, which_frame);
01232 int i;
01233 if (!ts) {
01234 for (i=0; i<num; i++) 
01235 if (flgs[i]) data[i] = 0.0;
01236 return 0; 
01237 }
01238 for (i=0; i<num; i++) { 
01239 if (flgs[i]) { 
01240 data[i] = ts->pos[3L*i + 1L]; 
01241 } 
01242 } 
01243 return 1; 
01244 }
01245 
01246 static int atomsel_zpos(void *v, int num, double *data, int *flgs) {
01247 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01248 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01249 Timestep *ts = selframe(atom_sel_mol, which_frame);
01250 int i;
01251 if (!ts) {
01252 for (i=0; i<num; i++) 
01253 if (flgs[i]) data[i] = 0.0;
01254 return 0; 
01255 }
01256 for (i=0; i<num; i++) { 
01257 if (flgs[i]) { 
01258 data[i] = ts->pos[3L*i + 2L]; 
01259 } 
01260 } 
01261 return 1; 
01262 }
01263 
01264 static int atomsel_ufx(void *v, int num, double *data, int *flgs) {
01265 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01266 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01267 Timestep *ts = selframe(atom_sel_mol, which_frame);
01268 int i;
01269 if (!ts || !ts->force) {
01270 for (i=0; i<num; i++) 
01271 if (flgs[i]) data[i] = 0.0;
01272 return 0; 
01273 }
01274 for (i=0; i<num; i++) { 
01275 if (flgs[i]) { 
01276 data[i] = ts->force[3L*i ];
01277 } 
01278 } 
01279 return 1; 
01280 }
01281 
01282 static int atomsel_ufy(void *v, int num, double *data, int *flgs) {
01283 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01284 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01285 Timestep *ts = selframe(atom_sel_mol, which_frame);
01286 int i;
01287 if (!ts || !ts->force) {
01288 for (i=0; i<num; i++) 
01289 if (flgs[i]) data[i] = 0.0;
01290 return 0; 
01291 }
01292 for (i=0; i<num; i++) { 
01293 if (flgs[i]) { 
01294 data[i] = ts->force[3L*i + 1L];
01295 } 
01296 } 
01297 return 1; 
01298 }
01299 
01300 static int atomsel_ufz(void *v, int num, double *data, int *flgs) {
01301 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01302 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01303 Timestep *ts = selframe(atom_sel_mol, which_frame);
01304 int i;
01305 if (!ts || !ts->force) {
01306 for (i=0; i<num; i++) 
01307 if (flgs[i]) data[i] = 0.0;
01308 return 0; 
01309 }
01310 for (i=0; i<num; i++) { 
01311 if (flgs[i]) { 
01312 data[i] = ts->force[3L*i + 2L];
01313 } 
01314 } 
01315 return 1; 
01316 }
01317 
01318 #define atomsel_get_vel(name, offset) \
01319 static int atomsel_##name(void *v, int num, double *data, int *flgs) { \
01320 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
01321 int which_frame = ((atomsel_ctxt *)v)->which_frame; \
01322 Timestep *ts = selframe(atom_sel_mol, which_frame); \
01323 int i; \
01324 if (!ts || !ts->vel) { \
01325 for (i=0; i<num; i++) \
01326 if (flgs[i]) data[i] = 0.0; \
01327 return 0; \
01328 } \
01329 for (i=0; i<num; i++) { \
01330 if (flgs[i]) { \
01331 data[i] = ts->vel[3L*i + (offset)]; \
01332 } \
01333 } \
01334 return 1; \
01335 } 
01336 
01337 atomsel_get_vel(vx, 0)
01338 atomsel_get_vel(vy, 1)
01339 atomsel_get_vel(vz, 2)
01340 
01341 static int atomsel_phi(void *v, int num, double *data, int *flgs) {
01342 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01343 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01344 Timestep *ts = selframe(atom_sel_mol, which_frame);
01345 int i;
01346 if (!ts) {
01347 for (i=0; i<num; i++) data[i] = 0;
01348 return 0; 
01349 } 
01350 const float *r = ts->pos; 
01351 for (i=0; i<num; i++) {
01352 if (!flgs[i]) continue;
01353 MolAtom *atom = atom_sel_mol->atom(i);
01354 int resid = atom->uniq_resid;
01355 int N = atom_sel_mol->find_atom_in_residue("N", resid);
01356 int CA = atom_sel_mol->find_atom_in_residue("CA", resid);
01357 int C = atom_sel_mol->find_atom_in_residue("C", resid);
01358 
01359 // Find the index of the C atom from the previous residue by searching
01360 // the atoms bonded to N for an atom named "C". 
01361 if (N < 0) {
01362 data[i] = 0;
01363 continue;
01364 }
01365 MolAtom *atomN = atom_sel_mol->atom(N);
01366 int cprev = -1;
01367 for (int j=0; j<atomN->bonds; j++) {
01368 int aindex = atomN->bondTo[j];
01369 int nameindex = atom_sel_mol->atom(aindex)->nameindex;
01370 if (!strcmp(atom_sel_mol->atomNames.name(nameindex), "C")) {
01371 cprev = aindex;
01372 break;
01373 }
01374 }
01375 if (cprev >= 0 && CA >= 0 && C >= 0) 
01376 data[i] = dihedral(r+3L*cprev, r+3L*N, r+3L*CA, r+3L*C);
01377 else
01378 data[i] = 0.0; 
01379 }
01380 return 0;
01381 }
01382 
01383 static int atomsel_psi(void *v, int num, double *data, int *flgs) {
01384 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01385 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01386 Timestep *ts = selframe(atom_sel_mol, which_frame);
01387 int i;
01388 if (!ts) {
01389 for (i=0; i<num; i++) data[i] = 0;
01390 return 0; 
01391 } 
01392 const float *r = ts->pos; 
01393 for (i=0; i<num; i++) {
01394 if (!flgs[i]) continue;
01395 MolAtom *atom = atom_sel_mol->atom(i);
01396 int resid = atom->uniq_resid;
01397 int N = atom_sel_mol->find_atom_in_residue("N", resid);
01398 int CA = atom_sel_mol->find_atom_in_residue("CA", resid);
01399 int C = atom_sel_mol->find_atom_in_residue("C", resid);
01400 
01401 // Find the index of the N atom from the next residue by searching the
01402 // atoms bonded to C for an atom named "N".
01403 if (C < 0) {
01404 data[i] = 0;
01405 continue;
01406 }
01407 MolAtom *atomC = atom_sel_mol->atom(C);
01408 int nextn = -1;
01409 for (int j=0; j<atomC->bonds; j++) {
01410 int aindex = atomC->bondTo[j];
01411 int nameindex = atom_sel_mol->atom(aindex)->nameindex;
01412 if (!strcmp(atom_sel_mol->atomNames.name(nameindex), "N")) {
01413 nextn = aindex;
01414 break;
01415 }
01416 }
01417 if (nextn >= 0 && N >= 0 && CA >= 0) 
01418 data[i] = dihedral(r+3L*N, r+3L*CA, r+3L*C, r+3L*nextn);
01419 else
01420 data[i] = 0.0; 
01421 }
01422 return 0;
01423 }
01424 
01425 static int atomsel_set_phi(void *v, int num, double *data, int *flgs) {
01426 // We rotate about the N-CA bond
01427 // 0. Get the current value of phi
01428 // 1. Translate, putting CA at the origin
01429 // 2. Compute the axis along N-CA
01430 // 3. Rotate just the N-terminus about the given axis
01431 // 4. Undo the translation
01432 
01433 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01434 SymbolTable *table = ((atomsel_ctxt *)v)->table;
01435 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01436 Timestep *ts = selframe(atom_sel_mol, which_frame);
01437 int i;
01438 if (!ts) {
01439 return 0; 
01440 } 
01441 float *r = ts->pos; 
01442 
01443 // Go through the residues; if the residue contains all the necessary atoms,
01444 // check to see if the CA of the residue is on. If it is, proceed with the
01445 // rotation. 
01446 for (i=0; i<atom_sel_mol->fragList.num(); i++) {
01447 Fragment *frag = atom_sel_mol->fragList[i];
01448 int nfragres = frag->residues.num();
01449 if (nfragres < 2) continue;
01450 int C = atom_sel_mol->find_atom_in_residue("C", frag->residues[0]);
01451 // Start at second residue since I need the previous residue for phi
01452 for (int ires = 1; ires < nfragres; ires++) {
01453 int resid = frag->residues[ires];
01454 int cprev = C;
01455 int N = atom_sel_mol->find_atom_in_residue("N", resid);
01456 int CA = atom_sel_mol->find_atom_in_residue("CA", resid);
01457 C = atom_sel_mol->find_atom_in_residue("C", resid);
01458 if (cprev < 0 || N < 0 || CA < 0 || C < 0) continue;
01459 if (!flgs[CA]) continue;
01460 float CApos[3], Npos[3], axis[3];
01461 vec_copy(CApos, r+3L*CA);
01462 vec_copy(Npos, r+3L*N);
01463 vec_sub(axis, Npos, CApos);
01464 vec_normalize(axis); 
01465 double oldphi = dihedral(r+3L*cprev, Npos, CApos, r+3L*C);
01466 // Select the N-terminal part of the fragment, which includes all
01467 // residues up to the current one, plus N and NH of the curent one.
01468 // I'm just going to create a new atom selection for now.
01469 
01470 AtomSel *nterm = new AtomSel(atom_sel_mol->app,
01471 table, atom_sel_mol->id());
01472 char buf[100];
01473 sprintf(buf, 
01474 "(fragment %d and residue < %d) or (residue %d and name N NH CA)",
01475 i, resid, resid);
01476 if (nterm->change(buf, atom_sel_mol) == AtomSel::NO_PARSE) {
01477 msgErr << "set phi: internal atom selection error" << sendmsg;
01478 msgErr << buf << sendmsg; 
01479 delete nterm;
01480 continue;
01481 }
01482 Matrix4 mat;
01483 mat.transvec(axis[0], axis[1], axis[2]);
01484 mat.rot((float) (data[CA]-oldphi), 'x');
01485 mat.transvecinv(axis[0], axis[1], axis[2]); 
01486 
01487 for (int j=0; j<num; j++) {
01488 if (nterm->on[j]) {
01489 float *pos = r+3L*j;
01490 vec_sub(pos, pos, CApos);
01491 mat.multpoint3d(pos, pos);
01492 vec_add(pos, pos, CApos); 
01493 }
01494 }
01495 delete nterm;
01496 }
01497 }
01498 return 0;
01499 }
01500 
01501 static int atomsel_set_psi(void *v, int num, double *data, int *flgs) {
01502 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
01503 SymbolTable *table = ((atomsel_ctxt *)v)->table;
01504 int which_frame = ((atomsel_ctxt *)v)->which_frame;
01505 Timestep *ts = selframe(atom_sel_mol, which_frame);
01506 if (!ts) {
01507 return 0; 
01508 }
01509 float *r = ts->pos;
01510 
01511 for (int i=0; i<atom_sel_mol->fragList.num(); i++) {
01512 Fragment *frag = atom_sel_mol->fragList[i];
01513 int nfragres = frag->residues.num();
01514 if (nfragres < 2) continue;
01515 int N = atom_sel_mol->find_atom_in_residue("N", frag->residues[nfragres-1]);
01516 for (int ires = nfragres-2; ires >= 0; ires--) {
01517 int resid = frag->residues[ires];
01518 int nextn = N;
01519 N = atom_sel_mol->find_atom_in_residue("N", resid);
01520 int CA = atom_sel_mol->find_atom_in_residue("CA", resid);
01521 int C = atom_sel_mol->find_atom_in_residue("C", resid);
01522 if (nextn < 0 || N < 0 || CA < 0 || C < 0) continue;
01523 if (!flgs[CA]) continue;
01524 float CApos[3], Cpos[3], axis[3];
01525 vec_copy(CApos, r+3L*CA);
01526 vec_copy(Cpos, r+3L*C);
01527 vec_sub(axis, Cpos, CApos);
01528 vec_normalize(axis);
01529 double oldpsi = dihedral(r+3L*N, CApos, Cpos, r+3L*nextn);
01530 
01531 AtomSel *cterm = new AtomSel(atom_sel_mol->app,
01532 table, atom_sel_mol->id());
01533 char buf[100];
01534 sprintf(buf,
01535 "(fragment %d and residue > %d) or (residue %d and name CA C O)",
01536 i, resid, resid);
01537 if (cterm->change(buf, atom_sel_mol) == AtomSel::NO_PARSE) {
01538 msgErr << "set psi: internal atom selection error" << sendmsg;
01539 msgErr << buf << sendmsg;
01540 delete cterm;
01541 continue;
01542 }
01543 
01544 Matrix4 mat;
01545 mat.transvec(axis[0], axis[1], axis[2]);
01546 mat.rot((float) (data[CA]-oldpsi), 'x');
01547 mat.transvecinv(axis[0], axis[1], axis[2]);
01548 
01549 for (int j=0; j<num; j++) {
01550 if (cterm->on[j]) {
01551 float *pos = r+3L*j;
01552 vec_sub(pos, pos, CApos);
01553 mat.multpoint3d(pos, pos);
01554 vec_add(pos, pos, CApos);
01555 }
01556 }
01557 delete cterm;
01558 }
01559 }
01560 return 0;
01561 }
01562 
01563 #define set_position_data(fctn, offset) \
01564 static int fctn(void *v, int num, double *data, int *flgs) \
01565 { \
01566 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
01567 int which_frame = ((atomsel_ctxt *)v)->which_frame; \
01568 Timestep *ts = selframe(atom_sel_mol, which_frame); \
01569 if (!ts) return 0; \
01570 for (int i=num-1; i>=0; i--) { \
01571 if (flgs[i]) { \
01572 ts->pos[3L*i + offset] = (float) data[i]; \
01573 } \
01574 } \
01575 return 1; \
01576 }
01577 
01578 
01579 set_position_data(atomsel_set_xpos, 0)
01580 set_position_data(atomsel_set_ypos, 1)
01581 set_position_data(atomsel_set_zpos, 2)
01582 
01583 #define set_force_data(fctn, offset) \
01584 static int fctn(void *v, int num, double *data, int *flgs) \
01585 { \
01586 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
01587 int which_frame = ((atomsel_ctxt *)v)->which_frame; \
01588 Timestep *ts = selframe(atom_sel_mol, which_frame); \
01589 if (!ts) return 0; \
01590 if (!ts->force) { \
01591 ts->force = new float[3L*num]; \
01592 memset(ts->force, 0, 3L*num*sizeof(float)); \
01593 } \
01594 for (int i=num-1; i>=0; i--) { \
01595 if (flgs[i]) { \
01596 ts->force[3L*i + offset] = (float) data[i]; \
01597 } \
01598 } \
01599 return 1; \
01600 }
01601 
01602 set_force_data(atomsel_set_ufx, 0)
01603 set_force_data(atomsel_set_ufy, 1)
01604 set_force_data(atomsel_set_ufz, 2)
01605 
01606 #define set_vel_data(fctn, offset) \
01607 static int fctn(void *v, int num, double *data, int *flgs) \
01608 { \
01609 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; \
01610 int which_frame = ((atomsel_ctxt *)v)->which_frame; \
01611 Timestep *ts = selframe(atom_sel_mol, which_frame); \
01612 if (!ts) return 0; \
01613 if (!ts->vel) { \
01614 ts->vel= new float[3L*num]; \
01615 memset(ts->vel, 0, 3L*num*sizeof(float)); \
01616 } \
01617 for (int i=num-1; i>=0; i--) { \
01618 if (flgs[i]) { \
01619 ts->vel[3L*i + offset] = (float) data[i]; \
01620 } \
01621 } \
01622 return 1; \
01623 }
01624 
01625 set_vel_data(atomsel_set_vx, 0)
01626 set_vel_data(atomsel_set_vy, 1)
01627 set_vel_data(atomsel_set_vz, 2)
01628 
01629 extern "C" {
01630 double atomsel_square(double x) { return x*x; }
01631 }
01632 
01633 // this is different than the previous. It allows me to search for a
01634 // given regex sequence. For instace, given the protein sequence
01635 // WAPDTYLVAPDAQD
01636 // the selection: sequence APDT
01637 // will select only the 2nd through 5th terms
01638 // the selection: sequence APD
01639 // will select 2nd through 4th, and 9th to 11th.
01640 // the selection: sequence "A.D"
01641 // will get 2-4 and 9-14
01642 // and so on.
01643 // If a residue name is not normal, it becomes an X
01644 
01645 // I am handed a list of strings for this selection
01646 // (eg, sequence APD "A.D" A to D)
01647 // Since there are no non-standard residue names (ie, no '*')
01648 // I'll interpret everything as a regex. Also, phrases like
01649 // "X to Y" are interpreted as "X.*Y"
01650 
01651 static int atomsel_sequence(void *v, int argc, const char **argv, int *types,
01652 int num, int *flgs)
01653 {
01654 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01655 int i;
01656 // make a temporary array for marking the selected atoms
01657 int *selected = new int[num];
01658 for (i=0; i<num; i++) {
01659 selected[i] = FALSE;
01660 }
01661 // make the list of regex'es
01662 JRegex **regex = (JRegex **) malloc( argc * sizeof(JRegex *));
01663 int num_r = 0;
01664 {
01665 JString pattern;
01666 for (i=0; i<argc; i++) {
01667 pattern = argv[i];
01668 if (types[i] >= 3) { // get the next term (if a "to" element)
01669 pattern += ".*";
01670 pattern += argv[++i];
01671 }
01672 regex[num_r] = new JRegex(pattern);
01673 num_r ++;
01674 }
01675 } // constructed the regex array
01676 if (num_r == 0) {
01677 return 0;
01678 }
01679 
01680 // construct a list of sequences from each protein (pfraglist)
01681 // and nucleic acid (nfragList)
01682 for (int fragcount=0; fragcount <2; fragcount++) {
01683 int pcount = atom_sel_mol->pfragList.num();
01684 int ncount = atom_sel_mol->nfragList.num();
01685 for (i=0; i< (fragcount == 0 ? pcount : ncount); i++) {
01686 int size = (fragcount == 0 ? atom_sel_mol->pfragList[i]->num()
01687 : atom_sel_mol->nfragList[i]->num());
01688 char *s = new char[size+1]; // so that it can be NULL-terminated
01689 char *t = s;
01690 int *mark = new int[size];
01691 
01692 for (int j=0; j<size; j++) {
01693 int residuenum = ((fragcount == 0) ?
01694 (*atom_sel_mol->pfragList[i])[j] : 
01695 (*atom_sel_mol->nfragList[i])[j]);
01696 int atomnum = (atom_sel_mol->residueList[residuenum]->atoms[0]);
01697 MolAtom *atom = atom_sel_mol->atom(atomnum);
01698 const char *resname = (atom_sel_mol->resNames).name(atom->resnameindex);
01699 mark[j] = FALSE;
01700 if (fragcount == 0) {
01701 // protein translations
01702 // PDB names
01703 if (!strcasecmp( resname, "GLY")) {*t++ = 'G'; continue;}
01704 if (!strcasecmp( resname, "ALA")) {*t++ = 'A'; continue;}
01705 if (!strcasecmp( resname, "VAL")) {*t++ = 'V'; continue;}
01706 if (!strcasecmp( resname, "PHE")) {*t++ = 'F'; continue;}
01707 if (!strcasecmp( resname, "PRO")) {*t++ = 'P'; continue;}
01708 if (!strcasecmp( resname, "MET")) {*t++ = 'M'; continue;}
01709 if (!strcasecmp( resname, "ILE")) {*t++ = 'I'; continue;}
01710 if (!strcasecmp( resname, "LEU")) {*t++ = 'L'; continue;}
01711 if (!strcasecmp( resname, "ASP")) {*t++ = 'D'; continue;}
01712 if (!strcasecmp( resname, "GLU")) {*t++ = 'E'; continue;}
01713 if (!strcasecmp( resname, "LYS")) {*t++ = 'K'; continue;}
01714 if (!strcasecmp( resname, "ARG")) {*t++ = 'R'; continue;}
01715 if (!strcasecmp( resname, "SER")) {*t++ = 'S'; continue;}
01716 if (!strcasecmp( resname, "THR")) {*t++ = 'T'; continue;}
01717 if (!strcasecmp( resname, "TYR")) {*t++ = 'Y'; continue;}
01718 if (!strcasecmp( resname, "HIS")) {*t++ = 'H'; continue;}
01719 if (!strcasecmp( resname, "CYS")) {*t++ = 'C'; continue;}
01720 if (!strcasecmp( resname, "ASN")) {*t++ = 'N'; continue;}
01721 if (!strcasecmp( resname, "GLN")) {*t++ = 'Q'; continue;}
01722 if (!strcasecmp( resname, "TRP")) {*t++ = 'W'; continue;}
01723 // CHARMM names
01724 if (!strcasecmp( resname, "HSE")) {*t++ = 'H'; continue;}
01725 if (!strcasecmp( resname, "HSD")) {*t++ = 'H'; continue;}
01726 if (!strcasecmp( resname, "HSP")) {*t++ = 'H'; continue;}
01727 // AMBER names
01728 if (!strcasecmp( resname, "CYX")) {*t++ = 'C'; continue;}
01729 } else {
01730 // nucleic acid translations
01731 if (!strcasecmp( resname, "ADE")) {*t++ = 'A'; continue;}
01732 if (!strcasecmp( resname, "A")) {*t++ = 'A'; continue;}
01733 if (!strcasecmp( resname, "THY")) {*t++ = 'T'; continue;}
01734 if (!strcasecmp( resname, "T")) {*t++ = 'T'; continue;}
01735 if (!strcasecmp( resname, "CYT")) {*t++ = 'C'; continue;}
01736 if (!strcasecmp( resname, "C")) {*t++ = 'C'; continue;}
01737 if (!strcasecmp( resname, "GUA")) {*t++ = 'G'; continue;}
01738 if (!strcasecmp( resname, "G")) {*t++ = 'G'; continue;}
01739 if (!strcasecmp( resname, "URA")) {*t++ = 'U'; continue;}
01740 if (!strcasecmp( resname, "U")) {*t++ = 'U'; continue;}
01741 }
01742 // then I have no idea
01743 *t++ = 'X';
01744 } // end loop 'j'; constructed the sequence for this protein
01745 *t = 0; // terminate the string
01746 
01747 // msgInfo << "sequence " << i << " is: " << s << sendmsg;
01748 // which of the residues match the regex(es)?
01749 for (int r=0; r<num_r; r++) {
01750 int len, start = 0, offset;
01751 while ((offset = (regex[r]->search(s, strlen(s), len, start)))
01752 != -1) {
01753 // then there was a match from offset to offset+len
01754 // msgInfo << "start " << start << " offset " << offset << " len "
01755 // << len << sendmsg;
01756 for (int loop=offset; loop<offset+len; loop++) {
01757 mark[loop] = 1;
01758 }
01759 start = offset+len;
01760 }
01761 }
01762 
01763 // the list of selected residues is in mark
01764 // turn on the right atoms
01765 for (int marked=0; marked<size; marked++) {
01766 if (mark[marked]) {
01767 int residuenum = (fragcount == 0 ?
01768 (*atom_sel_mol->pfragList[i])[marked] :
01769 (*atom_sel_mol->nfragList[i])[marked]);
01770 for (int atomloop=0;
01771 atomloop< atom_sel_mol->residueList[residuenum]->
01772 atoms.num(); atomloop++) {
01773 selected[atom_sel_mol->residueList[residuenum]->
01774 atoms[atomloop]]=TRUE;
01775 }
01776 }
01777 }
01778 delete [] mark;
01779 delete [] s;
01780 } // end loop i over the fragments
01781 } // end loop 'fragcount'
01782 
01783 
01784 // get rid of the compiled regex's
01785 for (i=0; i<num_r; i++) {
01786 delete regex[i];
01787 }
01788 // copy the 'selected' array into 'flgs'
01789 for (i=0; i<num; i++) {
01790 flgs[i] = flgs[i] && selected[i];
01791 }
01792 return 1;
01793 }
01794 
01795 /************ support for RasMol selections ******************/
01797 // the full form of a primitive is (seems to be)
01798 // {[<resname>]}{<resid>}{:<chain>}{.<atom name>}
01799 // if resname is only alpha, the [] can be dropped
01800 // if chain is alpha, the : can be dropped
01801 // resname only contains * if it is the first one ?
01802 // ? cannot go in the resid
01803 // * can only replace the whole field
01804 static int atomsel_rasmol_primitive(void *v, int argc, const char **argv, int *,
01805 int num, int *flgs)
01806 {
01807 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
01808 SymbolTable *table = ((atomsel_ctxt *)v)->table;
01809 
01810 // for each word, (ignoring the quote flags)
01811 for (int word=0; word<argc; word++) {
01812 const char *rnm0 = argv[word]; // resname start
01813 const char *rnm1; // and end position
01814 const char *rid0, *rid1;
01815 if (*rnm0 == '*') {
01816 rnm1 = rnm0 + 1;
01817 rid0 = rnm1;
01818 } else if (*rnm0 == '[') {
01819 rnm0++;
01820 rnm1 = rnm0;
01821 while (*rnm1 && *rnm1 != ']') { // find trailing bracket
01822 rnm1 ++;
01823 }
01824 if (rnm1 == rnm0) { // for cases like [] and "["
01825 rid0 = rnm1;
01826 } else {
01827 if (*rnm1==']') { // for cases like [so4]
01828 rid0 = rnm1+1;
01829 } else { // for (incorrect) cases like [so4
01830 rid0 = rnm1;
01831 }
01832 }
01833 } else { // then must be alpha or ?
01834 rnm1 = rnm0;
01835 while (isalpha(*rnm1) || *rnm1 == '?') { // find first non-alpha
01836 rnm1++;
01837 }
01838 rid0 = rnm1;
01839 }
01840 // got the resname
01841 
01842 // parse the resid
01843 rid1 = rid0;
01844 if (*rid1 == '*') {
01845 rid1++;
01846 } else {
01847 while (isdigit(*rid1)) {
01848 rid1++;
01849 }
01850 }
01851 
01852 // if this is the : delimiter, skip over it
01853 const char *chn0, *chn1;
01854 if (*rid1 == ':') {
01855 chn0 = rid1 + 1;
01856 } else {
01857 chn0 = rid1;
01858 }
01859 
01860 // get the chain
01861 // seek the . or end of string
01862 chn1 = chn0;
01863 while (*chn1 && *chn1 != '.') {
01864 chn1++;
01865 }
01866 
01867 const char *nm0, *nm1;
01868 if (*chn1 == '.') {
01869 nm0 = chn1 + 1;
01870 } else {
01871 nm0 = chn1;
01872 }
01873 nm1 = nm0;
01874 // seek the end of string
01875 while (*nm1) {
01876 nm1++;
01877 }
01878 
01879 
01880 // save the info into strings
01881 JString resname, resid, chain, name;
01882 const char *s;
01883 for (s=rnm0; s<rnm1; s++) {
01884 resname += *s;
01885 }
01886 for (s=rid0; s<rid1; s++) {
01887 resid += *s;
01888 }
01889 for (s=chn0; s<chn1; s++) {
01890 chain += *s;
01891 }
01892 for (s=nm0; s<nm1; s++) {
01893 name += *s;
01894 }
01895 // msgInfo << "resname: " << (const char *) resname << sendmsg;
01896 // msgInfo << "resid: " << (const char *) resid << sendmsg;
01897 // msgInfo << "chain: " << (const char *) chain << sendmsg;
01898 // msgInfo << "name: " << (const char *) name << sendmsg;
01899 
01900 // convert to the VMD regex ( ? => .? and * => .*)
01901 // (however, if there is a * for the whole field, delete the field)
01902 if (resname == "*") resname = "";
01903 if (resid == "*") resid = "";
01904 if (chain == "*") chain = "";
01905 if (name == "*") name = "";
01906 resname.gsub("?", ".?"); resname.gsub("*", ".*");
01907 resid.gsub("?", ".?"); resid.gsub("*", ".*");
01908 chain.gsub("?", ".?"); chain.gsub("*", ".*");
01909 name.gsub("?", ".?"); name.gsub("*", ".*");
01910 // make everything upcase
01911 resname.upcase();
01912 resid.upcase();
01913 chain.upcase();
01914 name.upcase();
01915 
01916 // construct a new search
01917 JString search;
01918 if (resname != "") {
01919 search = "resname ";
01920 search += '"';
01921 search += resname;
01922 search += '"';
01923 }
01924 if (resid != "") {
01925 if (search != "") {
01926 search += " and resid ";
01927 } else {
01928 search = "resid ";
01929 }
01930 search += '"';
01931 search += resid;
01932 search += '"';
01933 }
01934 // if the chain length > 1, it is a segname
01935 int is_segname = chain.length() > 1;
01936 if (chain != "") {
01937 if (search != "") {
01938 search += (is_segname ? " and segname " : " and chain ");
01939 } else {
01940 search = (is_segname ? "segname " : "chain ");
01941 }
01942 search += '"';
01943 search += chain;
01944 search += '"';
01945 }
01946 if (name != "") {
01947 if (search != "") {
01948 search += " and name ";
01949 } else {
01950 search = "name ";
01951 }
01952 search += '"';
01953 search += name;
01954 search += '"';
01955 }
01956 msgInfo << "Search = " << search << sendmsg;
01957 
01958 if (search == "") {
01959 search = "all";
01960 }
01961 // and do the search
01962 AtomSel *atomSel = new AtomSel(atom_sel_mol->app,
01963 table, atom_sel_mol->id());
01964 if (atomSel->change((const char *) search, atom_sel_mol) ==
01965 AtomSel::NO_PARSE) {
01966 msgErr << "rasmol: cannot understand: " << argv[word] << sendmsg;
01967 delete atomSel;
01968 continue;
01969 }
01970 
01971 // save the results
01972 {
01973 for (int i=0; i<num; i++) {
01974 flgs[i] = flgs[i] && atomSel->on[i];
01975 }
01976 }
01977 delete atomSel;
01978 }
01979 return 1;
01980 }
01981 
01982 int atomsel_custom_singleword(void *v, int num, int *flgs) {
01983 SymbolTable *table = ((atomsel_ctxt *)v)->table;
01984 const char *singleword = ((atomsel_ctxt *)v)->singleword;
01985 if (!singleword) {
01986 msgErr << "Internal error, atom selection macro called without context"
01987 << sendmsg;
01988 return 0;
01989 }
01990 const char *macro = table->get_custom_singleword(singleword);
01991 if (!macro) {
01992 msgErr << "Internal error, atom selection macro has lost its defintion."
01993 << sendmsg;
01994 return 0;
01995 }
01996 // Create new AtomSel based on the macro
01997 DrawMolecule *mol = ((atomsel_ctxt *)v)->atom_sel_mol;
01998 AtomSel *atomSel = new AtomSel(mol->app, table, mol->id());
01999 atomSel->which_frame = ((atomsel_ctxt *)v)->which_frame;
02000 if (atomSel->change(macro, mol) == AtomSel::NO_PARSE) {
02001 msgErr << "Unable to parse macro: " << macro << sendmsg;
02002 delete atomSel;
02003 return 0;
02004 }
02005 for (int i=0; i<num; i++) {
02006 flgs[i] = flgs[i] && atomSel->on[i];
02007 }
02008 delete atomSel;
02009 return 1;
02010 }
02011 
02012 
02013 
02014 
02015 // These functions allow voxel values from volumetric data loaded in a molecule
02016 // to be used in atom selections. Three keywords are implemented: volN returns
02017 // the value of the voxel under the selected atom, interpvolN returns the 
02018 // voxel value interpolated from neighboring voxels, and volindexN returns the 
02019 // numerical index of the underlying voxel (I found the latter useful, however
02020 // its usefulness might be too marginal to be included into VMD?)
02021 
02022 // This is a hack because the volume names are hardcoded. Ideally, VMD should
02023 // invent an "array keyword" which would contain a parsable index (e.g. 
02024 // "vol#2"). Right now, the first M keywords are hard-coded, and any volN with
02025 // N greater than M will not work. This is because adding an array keyword
02026 // involves a considerable amount of work, but perhaps this should be considered
02027 // eventually...
02028 
02029 // NOTE: if a coordinate lies outside the volmap range, a value of NAN is
02030 // assigned. Since NAN != NAN, NAN atoms will automatically never be included 
02031 // in selections (which is what we want!). Furthermore, "$sel get vol0" would
02032 // return a Tcl-parsable NAN to easily allow the user to identify coords that
02033 // lie outside the valid range.
02034 
02035 #ifndef NAN //not a number
02036 const float NAN = sqrtf(-1.f); //need some kind of portable NAN definition
02037 #endif
02038 
02039 
02040 static int atomsel_volume_array(void *v, int num, double *data, int *flgs, int array_index) {
02041 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
02042 int which_frame = ((atomsel_ctxt *)v)->which_frame;
02043 Timestep *ts = selframe(atom_sel_mol, which_frame);
02044 const VolumetricData *voldata = atom_sel_mol->get_volume_data(array_index);
02045 int i;
02046 if (!ts || !voldata) {
02047 if (!ts) msgErr << "atomselect: non-existent timestep." << sendmsg;
02048 if (!voldata) msgErr << "atomselect: volumetric map volume#" << array_index << " does not exist." << sendmsg;
02049 for (i=0; i<num; i++) if (flgs[i]) data[i] = NAN;
02050 return 0; 
02051 }
02052 for (i=0; i<num; i++) { 
02053 if (flgs[i]) 
02054 data[i] = voldata->voxel_value_from_coord(ts->pos[3L*i], ts->pos[3L*i+1], ts->pos[3L*i+2]); 
02055 } 
02056 return 1; 
02057 }
02058 
02059 
02060 static int atomsel_interp_volume_array(void *v, int num, double *data, int *flgs, int array_index) {
02061 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
02062 int which_frame = ((atomsel_ctxt *)v)->which_frame;
02063 Timestep *ts = selframe(atom_sel_mol, which_frame);
02064 const VolumetricData *voldata = atom_sel_mol->get_volume_data(array_index);
02065 int i;
02066 if (!ts || !voldata) {
02067 if (!ts) msgErr << "atomselect: non-existent timestep." << sendmsg;
02068 if (!voldata) msgErr << "atomselect: volumetric map volume#" << array_index << " does not exist." << sendmsg;
02069 for (i=0; i<num; i++) if (flgs[i]) data[i] = NAN;
02070 return 0; 
02071 }
02072 for (i=0; i<num; i++) { 
02073 if (flgs[i]) 
02074 data[i] = voldata->voxel_value_interpolate_from_coord(ts->pos[3L*i], ts->pos[3L*i+1], ts->pos[3L*i+2]); 
02075 } 
02076 return 1; 
02077 }
02078 
02079 
02080 static int atomsel_gridindex_array(void *v, int num, double *data, int *flgs, int array_index) {
02081 DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol; 
02082 int which_frame = ((atomsel_ctxt *)v)->which_frame;
02083 Timestep *ts = selframe(atom_sel_mol, which_frame);
02084 const VolumetricData *voldata = atom_sel_mol->get_volume_data(array_index);
02085 int i;
02086 if (!ts || !voldata) {
02087 if (!ts) msgErr << "atomselect: non-existent timestep." << sendmsg;
02088 if (!voldata) msgErr << "atomselect: volumetric map volume#" << array_index << " does not exist." << sendmsg;
02089 for (i=0; i<num; i++) if (flgs[i]) data[i] = NAN;
02090 return 0; 
02091 }
02092 for (i=0; i<num; i++) { 
02093 if (flgs[i]) 
02094 data[i] = double(voldata->voxel_index_from_coord(ts->pos[3L*i], ts->pos[3L*i+1], ts->pos[3L*i+2]));
02095 } 
02096 return 1; 
02097 }
02098 
02099 
02100 static int atomsel_gridindex0(void *v, int num, double *data, int *flgs) { 
02101 return atomsel_gridindex_array(v, num, data, flgs, 0); 
02102 }
02103 static int atomsel_gridindex1(void *v, int num, double *data, int *flgs) { 
02104 return atomsel_gridindex_array(v, num, data, flgs, 1); 
02105 }
02106 static int atomsel_gridindex2(void *v, int num, double *data, int *flgs) { 
02107 return atomsel_gridindex_array(v, num, data, flgs, 2); 
02108 }
02109 static int atomsel_gridindex3(void *v, int num, double *data, int *flgs) { 
02110 return atomsel_gridindex_array(v, num, data, flgs, 3); 
02111 }
02112 static int atomsel_gridindex4(void *v, int num, double *data, int *flgs) { 
02113 return atomsel_gridindex_array(v, num, data, flgs, 4); 
02114 }
02115 static int atomsel_gridindex5(void *v, int num, double *data, int *flgs) { 
02116 return atomsel_gridindex_array(v, num, data, flgs, 5); 
02117 }
02118 static int atomsel_gridindex6(void *v, int num, double *data, int *flgs) { 
02119 return atomsel_gridindex_array(v, num, data, flgs, 6); 
02120 }
02121 static int atomsel_gridindex7(void *v, int num, double *data, int *flgs) { 
02122 return atomsel_gridindex_array(v, num, data, flgs, 7); 
02123 }
02124 
02125 
02126 // 'volume0'
02127 static int atomsel_volume0(void *v, int num, double *data, int *flgs) { 
02128 return atomsel_volume_array(v, num, data, flgs, 0); 
02129 }
02130 static int atomsel_volume1(void *v, int num, double *data, int *flgs) { 
02131 return atomsel_volume_array(v, num, data, flgs, 1); 
02132 }
02133 static int atomsel_volume2(void *v, int num, double *data, int *flgs) { 
02134 return atomsel_volume_array(v, num, data, flgs, 2); 
02135 }
02136 static int atomsel_volume3(void *v, int num, double *data, int *flgs) { 
02137 return atomsel_volume_array(v, num, data, flgs, 3); 
02138 }
02139 static int atomsel_volume4(void *v, int num, double *data, int *flgs) { 
02140 return atomsel_volume_array(v, num, data, flgs, 4); 
02141 }
02142 static int atomsel_volume5(void *v, int num, double *data, int *flgs) { 
02143 return atomsel_volume_array(v, num, data, flgs, 5); 
02144 }
02145 static int atomsel_volume6(void *v, int num, double *data, int *flgs) { 
02146 return atomsel_volume_array(v, num, data, flgs, 6); 
02147 }
02148 static int atomsel_volume7(void *v, int num, double *data, int *flgs) { 
02149 return atomsel_volume_array(v, num, data, flgs, 7); 
02150 }
02151 
02152 
02153 // 'interpolated_volume0'
02154 static int atomsel_interp_volume0(void *v, int num, double *data, int *flgs) { 
02155 return atomsel_interp_volume_array(v, num, data, flgs, 0); 
02156 }
02157 static int atomsel_interp_volume1(void *v, int num, double *data, int *flgs) { 
02158 return atomsel_interp_volume_array(v, num, data, flgs, 1); 
02159 }
02160 static int atomsel_interp_volume2(void *v, int num, double *data, int *flgs) { 
02161 return atomsel_interp_volume_array(v, num, data, flgs, 2); 
02162 }
02163 static int atomsel_interp_volume3(void *v, int num, double *data, int *flgs) { 
02164 return atomsel_interp_volume_array(v, num, data, flgs, 3); 
02165 }
02166 static int atomsel_interp_volume4(void *v, int num, double *data, int *flgs) { 
02167 return atomsel_interp_volume_array(v, num, data, flgs, 4); 
02168 }
02169 static int atomsel_interp_volume5(void *v, int num, double *data, int *flgs) { 
02170 return atomsel_interp_volume_array(v, num, data, flgs, 5); 
02171 }
02172 static int atomsel_interp_volume6(void *v, int num, double *data, int *flgs) { 
02173 return atomsel_interp_volume_array(v, num, data, flgs, 6); 
02174 }
02175 static int atomsel_interp_volume7(void *v, int num, double *data, int *flgs) { 
02176 return atomsel_interp_volume_array(v, num, data, flgs, 7); 
02177 }
02178 
02179 
02180 
02181 
02182 void atomSelParser_init(SymbolTable *atomSelParser) {
02183 atomSelParser->add_keyword("name", atomsel_name, atomsel_set_name);
02184 atomSelParser->add_keyword("type", atomsel_type, atomsel_set_type);
02185 
02186 // XXX probably only use these two for testing of BaseMolecule::analyze()
02187 atomSelParser->add_keyword("backbonetype", atomsel_backbonetype, NULL);
02188 atomSelParser->add_keyword("residuetype", atomsel_residuetype, NULL);
02189 
02190 atomSelParser->add_keyword("index", atomsel_index, NULL);
02191 atomSelParser->add_keyword("serial", atomsel_serial, NULL);
02192 atomSelParser->add_keyword("atomicnumber", atomsel_atomicnumber, atomsel_set_atomicnumber);
02193 atomSelParser->add_keyword("element", atomsel_element, atomsel_set_element);
02194 atomSelParser->add_keyword("residue", atomsel_residue, NULL);
02195 atomSelParser->add_keyword("resname", atomsel_resname, atomsel_set_resname);
02196 atomSelParser->add_keyword("altloc", atomsel_altloc, atomsel_set_altloc);
02197 atomSelParser->add_keyword("resid", atomsel_resid, atomsel_set_resid);
02198 atomSelParser->add_keyword("insertion", atomsel_insertion, NULL);
02199 atomSelParser->add_keyword("chain", atomsel_chain, atomsel_set_chain);
02200 atomSelParser->add_keyword("segname", atomsel_segname, atomsel_set_segname);
02201 atomSelParser->add_keyword("segid", atomsel_segname, atomsel_set_segname);
02202 
02203 atomSelParser->add_singleword("all", atomsel_all, NULL);
02204 atomSelParser->add_singleword("none", atomsel_none, NULL);
02205 
02206 atomSelParser->add_keyword("fragment", atomsel_fragment, NULL);
02207 atomSelParser->add_keyword("pfrag", atomsel_pfrag, NULL);
02208 atomSelParser->add_keyword("nfrag", atomsel_nfrag, NULL);
02209 atomSelParser->add_keyword("numbonds", atomsel_numbonds, NULL);
02210 
02211 atomSelParser->add_singleword("backbone", atomsel_backbone, NULL);
02212 atomSelParser->add_singleword("sidechain", 
02213 atomsel_sidechain, NULL);
02214 atomSelParser->add_singleword("protein", atomsel_protein, NULL);
02215 atomSelParser->add_singleword("nucleic", atomsel_nucleic, NULL);
02216 atomSelParser->add_singleword("water", atomsel_water, NULL);
02217 atomSelParser->add_singleword("waters", atomsel_water, NULL);
02218 atomSelParser->add_singleword("vmd_fast_hydrogen", atomsel_hydrogen, NULL);
02219 
02220 // secondary structure functions
02221 atomSelParser->add_singleword("helix", 
02222 atomsel_helix, atomsel_set_helix);
02223 atomSelParser->add_singleword("alpha_helix", 
02224 atomsel_alpha_helix, atomsel_set_helix);
02225 atomSelParser->add_singleword("helix_3_10", 
02226 atomsel_3_10_helix, atomsel_set_3_10_helix);
02227 atomSelParser->add_singleword("pi_helix", 
02228 atomsel_pi_helix, atomsel_set_pi_helix);
02229 atomSelParser->add_singleword("sheet", atomsel_sheet, atomsel_set_sheet);
02230 atomSelParser->add_singleword("betasheet", atomsel_sheet, atomsel_set_sheet);
02231 atomSelParser->add_singleword("beta_sheet",atomsel_sheet, atomsel_set_sheet);
02232 atomSelParser->add_singleword("extended_beta", atomsel_extended_sheet, 
02233 atomsel_set_sheet);
02234 atomSelParser->add_singleword("bridge_beta", atomsel_bridge_sheet, 
02235 atomsel_set_bridge_sheet);
02236 atomSelParser->add_singleword("turn", atomsel_turn, atomsel_set_turn);
02237 atomSelParser->add_singleword("coil", atomsel_coil, atomsel_set_coil);
02238 atomSelParser->add_keyword("structure",atomsel_sstruct, atomsel_set_sstruct);
02239 
02240 #if defined(VMDWITHCARBS)
02241 atomSelParser->add_keyword("pucker", atomsel_pucker, NULL);
02242 #endif
02243 
02244 atomSelParser->add_keyword("user", atomsel_user, atomsel_set_user);
02245 atomSelParser->add_keyword("user2", atomsel_user2, atomsel_set_user2);
02246 atomSelParser->add_keyword("user3", atomsel_user3, atomsel_set_user3);
02247 atomSelParser->add_keyword("user4", atomsel_user4, atomsel_set_user4);
02248 
02249 atomSelParser->add_keyword("x", atomsel_xpos, atomsel_set_xpos);
02250 atomSelParser->add_keyword("y", atomsel_ypos, atomsel_set_ypos);
02251 atomSelParser->add_keyword("z", atomsel_zpos, atomsel_set_zpos);
02252 atomSelParser->add_keyword("vx", atomsel_vx, atomsel_set_vx);
02253 atomSelParser->add_keyword("vy", atomsel_vy, atomsel_set_vy);
02254 atomSelParser->add_keyword("vz", atomsel_vz, atomsel_set_vz);
02255 atomSelParser->add_keyword("ufx", atomsel_ufx, atomsel_set_ufx);
02256 atomSelParser->add_keyword("ufy", atomsel_ufy, atomsel_set_ufy);
02257 atomSelParser->add_keyword("ufz", atomsel_ufz, atomsel_set_ufz);
02258 atomSelParser->add_keyword("phi", atomsel_phi, atomsel_set_phi);
02259 atomSelParser->add_keyword("psi", atomsel_psi, atomsel_set_psi);
02260 atomSelParser->add_keyword("radius", atomsel_radius, 
02261 atomsel_set_radius);
02262 atomSelParser->add_keyword("mass", atomsel_mass, atomsel_set_mass);
02263 atomSelParser->add_keyword("charge", atomsel_charge, 
02264 atomsel_set_charge);
02265 atomSelParser->add_keyword("beta", atomsel_beta, atomsel_set_beta);
02266 atomSelParser->add_keyword("occupancy", 
02267 atomsel_occupancy, atomsel_set_occupancy);
02268 
02269 atomSelParser->add_keyword("flag0", atomsel_flag0, atomsel_set_flag0);
02270 atomSelParser->add_keyword("flag1", atomsel_flag1, atomsel_set_flag1);
02271 atomSelParser->add_keyword("flag2", atomsel_flag2, atomsel_set_flag2);
02272 atomSelParser->add_keyword("flag3", atomsel_flag3, atomsel_set_flag3);
02273 atomSelParser->add_keyword("flag4", atomsel_flag4, atomsel_set_flag4);
02274 atomSelParser->add_keyword("flag5", atomsel_flag5, atomsel_set_flag5);
02275 atomSelParser->add_keyword("flag6", atomsel_flag6, atomsel_set_flag6);
02276 atomSelParser->add_keyword("flag7", atomsel_flag7, atomsel_set_flag7);
02277 
02278 atomSelParser->add_stringfctn("sequence", atomsel_sequence);
02279 atomSelParser->add_stringfctn("rasmol", 
02280 atomsel_rasmol_primitive);
02281 
02282 // three letters for resname, 1 or more letters for resid
02284 // atomSelParser->add_singleword("[a-zA-Z][a-zA-Z][a-zA-Z][0-9]+",
02285 // "resnameID", atomsel_resname_resid,
02286 // NULL);
02287 
02288 // and a few functions for good measure
02289 // Note: These functions must return the same output for given input.
02290 // Functions like rand() will break some optimizations in the 
02291 // atom selection code otherwise.
02292 atomSelParser->add_cfunction("sqr", atomsel_square);
02293 atomSelParser->add_cfunction("sqrt", sqrt);
02294 atomSelParser->add_cfunction("abs", fabs);
02295 atomSelParser->add_cfunction("floor", floor);
02296 atomSelParser->add_cfunction("ceil", ceil);
02297 atomSelParser->add_cfunction("sin", sin);
02298 atomSelParser->add_cfunction("cos", cos);
02299 atomSelParser->add_cfunction("tan", tan);
02300 atomSelParser->add_cfunction("atan", atan);
02301 atomSelParser->add_cfunction("asin", asin);
02302 atomSelParser->add_cfunction("acos", acos);
02303 atomSelParser->add_cfunction("sinh", sinh);
02304 atomSelParser->add_cfunction("cosh", cosh);
02305 atomSelParser->add_cfunction("tanh", tanh);
02306 atomSelParser->add_cfunction("exp", exp);
02307 atomSelParser->add_cfunction("log", log);
02308 atomSelParser->add_cfunction("log10", log10);
02309 
02310 
02311 
02312 atomSelParser->add_keyword("volindex0", atomsel_gridindex0, NULL);
02313 atomSelParser->add_keyword("volindex1", atomsel_gridindex1, NULL);
02314 atomSelParser->add_keyword("volindex2", atomsel_gridindex2, NULL);
02315 atomSelParser->add_keyword("volindex3", atomsel_gridindex3, NULL);
02316 atomSelParser->add_keyword("volindex4", atomsel_gridindex4, NULL);
02317 atomSelParser->add_keyword("volindex5", atomsel_gridindex5, NULL);
02318 atomSelParser->add_keyword("volindex6", atomsel_gridindex6, NULL);
02319 atomSelParser->add_keyword("volindex7", atomsel_gridindex7, NULL);
02320 atomSelParser->add_keyword("vol0", atomsel_volume0, NULL);
02321 atomSelParser->add_keyword("vol1", atomsel_volume1, NULL);
02322 atomSelParser->add_keyword("vol2", atomsel_volume2, NULL);
02323 atomSelParser->add_keyword("vol3", atomsel_volume3, NULL);
02324 atomSelParser->add_keyword("vol4", atomsel_volume4, NULL);
02325 atomSelParser->add_keyword("vol5", atomsel_volume5, NULL);
02326 atomSelParser->add_keyword("vol6", atomsel_volume6, NULL);
02327 atomSelParser->add_keyword("vol7", atomsel_volume7, NULL);
02328 atomSelParser->add_keyword("interpvol0", atomsel_interp_volume0, NULL);
02329 atomSelParser->add_keyword("interpvol1", atomsel_interp_volume1, NULL);
02330 atomSelParser->add_keyword("interpvol2", atomsel_interp_volume2, NULL);
02331 atomSelParser->add_keyword("interpvol3", atomsel_interp_volume3, NULL);
02332 atomSelParser->add_keyword("interpvol4", atomsel_interp_volume4, NULL);
02333 atomSelParser->add_keyword("interpvol5", atomsel_interp_volume5, NULL);
02334 atomSelParser->add_keyword("interpvol6", atomsel_interp_volume6, NULL);
02335 atomSelParser->add_keyword("interpvol7", atomsel_interp_volume7, NULL);
02336 }
02337 
02338 
02340 // constructor; parse string and see if OK
02341 AtomSel::AtomSel(VMDApp *vmdapp, SymbolTable *st, int id)
02342 : ID(id) {
02343 
02344 // initialize variables
02345 app = vmdapp;
02346 table = st;
02347 selected = NO_PARSE;
02348 firstsel = 0;
02349 lastsel = NO_PARSE;
02350 on = NULL;
02351 on256 = NULL;
02352 cmdStr = NULL;
02353 tree = NULL;
02354 num_atoms = 0;
02355 num_atoms256 = 0;
02356 which_frame = TS_NOW;
02357 do_update = 0;
02358 }
02359 
02360 // destructor; free up space and make all pointers invalid
02361 AtomSel::~AtomSel(void) {
02362 table = NULL;
02363 num_atoms = 0;
02364 num_atoms256 = 0;
02365 selected = NO_PARSE;
02366 firstsel = 0;
02367 lastsel = NO_PARSE;
02368 delete [] on;
02369 on = NULL;
02370 delete [] on256;
02371 on256 = NULL;
02372 delete tree;
02373 delete [] cmdStr;
02374 }
02375 
02376 int AtomSel::change(const char *newcmd, DrawMolecule *mol) {
02377 if (newcmd) {
02378 ParseTree *newtree = table->parse(newcmd);
02379 if (!newtree) {
02380 return NO_PARSE;
02381 }
02382 delete [] cmdStr;
02383 cmdStr = stringdup(newcmd);
02384 delete tree;
02385 tree = newtree;
02386 }
02387 
02388 // and evaluate
02389 atomsel_ctxt context(table, mol, which_frame, NULL);
02390 
02391 // resize flags array if necessary
02392 if (num_atoms != mol->nAtoms || (on == NULL)) {
02393 if (on) delete [] on;
02394 on = new int[mol->nAtoms];
02395 num_atoms = mol->nAtoms;
02396 }
02397 
02398 tree->use_context(&context);
02399 int rc = tree->evaluate(mol->nAtoms, on);
02400 
02401 // store parse return code before we postprocess
02402 int ret_code = rc ? PARSE_SUCCESS : NO_EVAL;
02403 
02404 // find the first selected atom, the last selected atom,
02405 // and the total number of selected atoms
02406 selected = 0;
02407 firstsel = 0; // ensure that if nothing is selected, that we
02408 lastsel = -1; // trim subsequent loops to do no iterations
02409 
02410 // use high performance vectorized select analysis routine
02411 if (analyze_selection_aligned_dispatch(app->cpucaps, num_atoms, on, &firstsel, &lastsel, &selected))
02412 return ret_code;
02413 
02414 #if 0
02415 // (re)allocate and populate selection group acceleration data structure
02416 num_atoms256 = num_atoms >> 8; // natoms/256
02417 if (on256) delete [] on256;
02418 on256 = new unsigned int[num_atoms256];
02419 
02420 int firstblk = firstsel >> 8;
02421 int lastblk = lastsel >> 8;
02422 for (int blk=firstblk; blk<=lastblk; blk++) {
02423 // loop and test only the atoms in block[blk]
02424 int firstatom = blk << 8;
02425 int lastatom = ((blk+1) << 8) - 1;
02426 int blkcnt = 0;
02427 int i=0;
02428 
02429 on256 = 0; // no atoms selected initially
02430 for (i=firstatom; i<=lastatom; i++) {
02431 if (on[i]) {
02432 // set selection flag non-zero indicating that something within
02433 // this 256-atom block is selected...
02434 on256[blk] = 1 << 24; // reserve lower 3 bytes for indices/counter
02435 
02436 break; // once we've found a selected atom, we can test next blk
02437 }
02438 }
02439 
02440 #if 0
02441 // XXX a next-gen scheme could break up the 32-bit integer into
02442 // four 8-bit quantities that indicate the first-selected 
02443 // index (relative to the start of the block), the last-selected
02444 // index (relative to the start of the block), and the total number
02445 // of selected atoms in the block. Since we work with 256-atom
02446 // blocks, each unsigned sub-byte can fully index/count the 
02447 // contents of the block, allowing block-level optimizations 
02448 // identical in design to what we have done previously at the
02449 // level of entire selection arrays.
02450 on[i] |= (i & 0xFF) << 8; // first selected in block 
02451 for (; i<=lastatom; i++) {
02452 if (on[i]) {
02453 blkcnt++; // count selected atoms...
02454 }
02455 }
02456 on[i] |= (i & 0xFF) << 16; // last selected in block 
02457 on[i] |= (blkcnt & 0xFF); // count of selected atoms in block
02458 #endif
02459 
02460 }
02461 #endif
02462 
02463 return ret_code;
02464 }
02465 
02466 
02467 // return the current coordinates (or NULL if error)
02468 float *AtomSel::coordinates(MoleculeList *mlist) const {
02469 Timestep *ts = timestep(mlist);
02470 if (ts) return ts->pos;
02471 return NULL;
02472 }
02473 
02474 
02475 // return the current timestep (or NULL if error)
02476 Timestep *AtomSel::timestep(MoleculeList *mlist) const {
02477 DrawMolecule *mymol = mlist->mol_from_id(molid()); 
02478 
02479 if (!mymol) {
02480 msgErr << "No molecule" << sendmsg;
02481 return NULL; // no molecules
02482 }
02483 
02484 switch (which_frame) {
02485 case TS_LAST: 
02486 return mymol->get_last_frame();
02487 
02488 case TS_NOW: 
02489 return mymol->current();
02490 
02491 default:
02492 // if past end of coords return last coord
02493 if (!mymol->get_frame(which_frame)) { 
02494 return mymol->get_last_frame(); 
02495 }
02496 }
02497 
02498 return mymol->get_frame(which_frame);
02499 }
02500 
02501 
02502 int AtomSel::get_frame_value(const char *s, int *val) {
02503 *val = 1;
02504 if (!strcmp(s, "last")) {
02505 *val = TS_LAST;
02506 }
02507 if (!strcmp(s, "first")) {
02508 *val = 0;
02509 }
02510 if (!strcmp(s, "now")) {
02511 *val = TS_NOW;
02512 }
02513 if (*val == 1) {
02514 int new_frame = atoi(s);
02515 *val = new_frame;
02516 if (new_frame < 0) {
02517 return -1;
02518 }
02519 }
02520 return 0;
02521 }