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

AtomSel.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the           
00004  *cr                        University of Illinois                       
00005  *cr                         All Rights Reserved                        
00006  *cr                                                                   
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
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 }

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