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

MeasureSurface.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: MeasureSurface.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.11 $       $Date: 2019/09/27 05:36:06 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code to find surface atoms
00019  ***************************************************************************/
00020 
00021 #include <math.h>
00022 #include <stdlib.h>
00023 #include <stdio.h>
00024 #include "Measure.h"
00025 
00026 class Vector3D;
00027 class IndexList;
00028 class AtomData;
00029 class MolData;
00030 class AtomList;
00031 class AtomGrid;
00032 class GridList;
00033 class AtomGridCell;
00034 
00035 // class Vector3D
00036 // Operations on 3-D double vectors
00037 //
00038 class Vector3D {
00039 public:
00040    const Vector3D operator+(const Vector3D w ) const {
00041       Vector3D v;
00042       v.x = x + w.x;
00043       v.y = y + w.y;
00044       v.z = z + w.z;
00045       return v;
00046    }
00047    
00048    const Vector3D operator-(const Vector3D w ) const {
00049       Vector3D v;
00050       v.x = x - w.x;
00051       v.y = y - w.y;
00052       v.z = z - w.z;
00053       return v;
00054    }
00055    
00056    const Vector3D operator*(double s) const {
00057       Vector3D v;
00058       v.x = s*x;
00059       v.y = s*y;
00060       v.z = s*z;
00061       return v;
00062    }
00063    
00064    double length() const { return sqrt(x*x + y*y + z*z); }
00065    double length2() const { return x*x + y*y + z*z; }
00066    double x,y,z;
00067 };
00068 
00069 Vector3D operator*(double s, Vector3D v) {
00070    v.x *= s;
00071    v.y *= s;
00072    v.z *= s;
00073    return v;
00074 }
00075 
00076 Vector3D operator/(Vector3D v, double s) {
00077    v.x /= s;
00078    v.y /= s;
00079    v.z /= s;
00080    return v;
00081 }
00082 
00083 // class IndexList
00084 // A grow-able list of integers, for holding indices of atoms
00085 //
00086 class IndexList {
00087 public:
00088    IndexList() {
00089       num=0;
00090       maxnum=16;
00091       list = new int[16];
00092    }
00093    
00094    ~IndexList() {
00095       if (list != 0)
00096          delete [] list;
00097    }
00098    
00099    void add(const int val) {
00100       // If we need more space, allocate a new block that is
00101       //  1.5 times larger and copy everything over
00102       if (num == maxnum) {
00103          maxnum = (int)(maxnum * 1.5 + 1);
00104          int* oldlist = list;
00105          list = new int[maxnum];
00106          int i;
00107          for(i=0;i<num;i++) {
00108             list[i] = oldlist[i];
00109          }
00110          delete [] oldlist;
00111       }
00112       
00113       // We should have enough space now, add the value
00114       list[num] = val;
00115       num++;
00116    }
00117    
00118    int size() {
00119       return num;
00120    }
00121    
00122    int get(const int i) const { return list[i]; }
00123    
00124 private:
00125    int num;
00126    int maxnum;
00127    int* list;
00128 };
00129 
00130 // class AtomData
00131 // Store the atom index and coordinate in one object
00132 //
00133 class AtomData {
00134 public:
00135    int index;
00136    Vector3D coord;
00137 };
00138 
00139 // class MolData
00140 // Class for storing and reading in the atom coordinates being processed.
00141 // It also contains the routine for determining surface atoms, once the
00142 // vacuum-grid has been built
00143 //
00144 class MolData {
00145 public:
00146    MolData() {
00147       coords = 0;
00148    }
00149    
00150    ~MolData() {
00151       if ( coords != 0) {
00152          delete [] coords;
00153       }
00154    }
00155    
00156    static MolData* readFile(const char* fname,
00157                             double gridspacing, double radius, double threshold);
00158 
00159    static MolData* readData(const int *indx, const float *coords, const int count,
00160                             const Vector3D o,
00161                             const Vector3D a, const Vector3D b, const Vector3D c,
00162                             const double gridspacing, const double radius, 
00163                             const double threshold,
00164                             const bool periodic);
00165    
00166    void wrapCoords(const AtomGrid* grid);
00167    
00168    AtomList* findAtomsGrid(const AtomGrid* grid);
00169    double getRadius() const { return radius; }
00170    double getActualRadius() const { return actual_radius; }
00171    
00172    double gridspacing;
00173    double radius;
00174    double threshold;
00175    double actual_radius;
00176    Vector3D origin;
00177    Vector3D basisa;
00178    Vector3D basisb;
00179    Vector3D basisc;
00180    int count;
00181    AtomData* coords;
00182    bool periodic;
00183 };
00184 
00185 // class AtomList
00186 // Stores the atoms that are part of the shell, and writes their indices out
00187 // at the end
00188 //
00189 class AtomList {
00190 public:
00191    AtomList(const int sz) {
00192       size = sz;
00193       index_list = new int[sz];
00194       count = 0;
00195    }
00196    
00197    ~AtomList() {
00198       if (index_list != 0) {
00199          delete [] index_list;
00200       }
00201    }
00202    
00203    void add(const int index) { index_list[count] = index; count++; }
00204    
00205    void storeFile(const char* fname) const 
00206    {
00207       FILE* outf = fopen(fname,"w");
00208       if (outf == NULL) {
00209          printf("Couldn't open file %s\n",fname);
00210          exit(-1);
00211       }
00212       
00213       int i;
00214       for(i=0; i < count; i++) {
00215          if (fprintf(outf,"%d\n",index_list[i]) < 0) {
00216             printf("Store file error %d %s\n",i,fname);
00217             exit(-2);
00218          }
00219       }
00220       
00221       fclose(outf);
00222       return;
00223    }
00224    
00225    void storeData(int **surf, int *surf_count) const 
00226    {
00227       *surf_count = count;
00228       *surf = new int[count];
00229       int i;
00230       for(i=0; i < count; i++) {
00231          (*surf)[i] = index_list[i];
00232       }
00233       
00234       return;
00235    }
00236    
00237 private:
00238    int* index_list;
00239    int count;
00240    int size;
00241 };
00242 
00243 // class GridList
00244 // Store a list of grid points, with each grid point identified as a 3-D vector
00245 // of ints
00246 //
00247 class GridList {
00248 public:
00249    GridList(AtomGrid* thegrid, const int gridn)
00250    {
00251       n = gridn;
00252       if (n != 0)
00253          items = new Vector3D[n];
00254       else
00255          items = 0;
00256       
00257       indx = 0;
00258       grid = thegrid;
00259       shuffle_vec = scramble_array(n);
00260    }
00261    
00262    ~GridList() {
00263       if (items != 0)
00264          delete [] items;
00265       if (shuffle_vec != 0)
00266          delete [] shuffle_vec;
00267    }
00268    
00269    int add(const int i, const int j, const int k);
00270    Vector3D get(const int i) const {
00271       if (i < 0 || i >= n)
00272          return Vector3D(); // XXX: this used to throw an (unhandled) exception
00273       
00274       return items[i];
00275    }
00276    
00277    int getNum() const { return indx; }
00278    
00279 private:
00280    int n;
00281    int indx;
00282    AtomGrid* grid;
00283    int* shuffle_vec;
00284    Vector3D* items;
00285    static int* scramble_array(const int n);
00286 };
00287 
00288 // class AtomGridCell
00289 // Divide the grid into a mesh of cubic cells. Each cell contains all the grid
00290 // points for that region, and can produce a list of "set" grid points on
00291 // demand
00292 //
00293 class AtomGridCell {
00294 public:
00295    AtomGridCell(AtomGrid* in_ag, int in_mina, int in_maxa, int in_minb, 
00296                 int in_maxb, int in_minc, int in_maxc)
00297    {
00298       //    std::cout << "New AtomGridCell " << mina << "," << maxa << ","
00299       //      << minb << "," << maxb << ","
00300       //      << minc << "," << maxc << ","
00301       //      << std::endl;
00302       
00303       ag = in_ag;
00304       mina = in_mina;
00305       maxa = in_maxa;
00306       minb = in_minb;
00307       maxb = in_maxb;
00308       minc = in_minc;
00309       maxc = in_maxc;
00310       
00311       na = maxa-mina+1;
00312       nb = maxb-minb+1;
00313       nc = maxc-minc+1;
00314       
00315       noneset = true;
00316       getlistcalled = false;
00317       cellmap = 0;
00318       count=0;
00319       gridlist=0;
00320    }
00321    
00322    ~AtomGridCell()
00323    {
00324       if (cellmap != 0)
00325          delete [] cellmap;
00326    }
00327    
00328    int set(const int i, const int j, const int k);
00329    bool get(const int i, const int j, const int k) const;
00330    const GridList* get_list();
00331    
00332 private:
00333    void build_list();
00334    AtomGrid* ag;
00335    int mina;
00336    int maxa;
00337    int minb;
00338    int maxb;
00339    int minc;
00340    int maxc;
00341    int na;
00342    int nb;
00343    int nc;
00344    bool noneset;
00345    bool getlistcalled;
00346    bool* cellmap;
00347    int count;
00348    GridList* gridlist;
00349 };
00350 
00351 // class AtomGrid
00352 // Build and maintain a grid of boolean values for the entire space
00353 //
00354 class AtomGrid {
00355 public:
00356    AtomGrid(const MolData* mdata);
00357    
00358    ~AtomGrid() {
00359       int i,j,k;
00360       for(i=0; i < cna; i++)
00361          for(j=0; j < cnb; j++)
00362             for(k=0; k < cnc; k++) {
00363                delete cellgrid[((i*cnb)+j)*cnc+k];
00364             }
00365       delete [] cellgrid;
00366       delete [] celli;
00367       delete [] cellj;
00368       delete [] cellk;
00369    }
00370    
00371    // get(const int i, const int j, const int k)
00372    // Returns the state (true/false) of grid point (i,j,k)
00373    //
00374    bool get(const int i, const int j, const int k) const 
00375    {
00376       // Find which cell we need, then get the value from that cell
00377       bool ret;
00378       if (i >= 0 && i < na && j >= 0 && j < nb && k >= 0 && k < nc) {
00379          const int ii = celli[i];
00380          const int jj = cellj[j];
00381          const int kk = cellk[k];
00382          const int indx = ((ii*cnb)+jj)*cnc+kk;
00383          ret = cellgrid[indx]->get(i,j,k);
00384       } else {
00385          ret = false;
00386       }
00387       return ret;
00388    }
00389    
00390    // set(const int i, const int j, const int k)
00391    // Set grid point (i,j,k) true. Default is false
00392    //
00393    int set(const int i, const int j, const int k) {
00394       //    std::cout << "set called " << i << "," << j << "." << k << std::endl;
00395       if ( i < 0 || i >= na || j < 0 || j >= nb || k < 0 || k >= nc )
00396          return -1;
00397       
00398       const int ii = celli[i];
00399       const int jj = cellj[j];
00400       const int kk = cellk[k];
00401       return cellgrid[((ii*cnb)+jj)*cnc+kk]->set(i,j,k);
00402    }
00403    
00404    // get_ijk(const Vector3D vv) const
00405    // Convert the real coordinates vv to grid coordinate space.
00406    // The grid coordinates may have a decimal component, indicating
00407    // where it falls within the grid block.
00408    //
00409    Vector3D get_ijk(const Vector3D vv) const 
00410    {
00411       Vector3D ijk;
00412       const Vector3D v = vv - origin;
00413       //    std::cout << "Pos = " << v.x << " " << v.y << " " << v.z
00414       //               << "---" << vv.x << " " << vv.y << " " << vv.z << std::endl;
00415       
00416       ijk.x = (inv_abc[0][0] * v.x + inv_abc[0][1] * v.y + inv_abc[0][2] * v.z)
00417          * na;
00418       ijk.y = (inv_abc[1][0] * v.x + inv_abc[1][1] * v.y + inv_abc[1][2] * v.z)
00419          * nb;
00420       ijk.z = (inv_abc[2][0] * v.x + inv_abc[2][1] * v.y + inv_abc[2][2] * v.z)
00421          * nc;
00422       
00423       return ijk;
00424    }
00425    
00426    // get_cijk(const Vector3D vv, int& i, int& j int& k) const
00427    // Convert the real coordinates vv to grid coordinate space, then determine
00428    // the coordinates of the cell which owns that grid point, and return
00429    // those coordinates in i, j, k
00430    //
00431    void get_cijk(const Vector3D vv, int& i, int& j, int& k) const 
00432    {
00433       Vector3D ijk = get_ijk(vv);
00434       
00435       ijk.x = fmod(ijk.x,na);
00436       if (ijk.x < 0)
00437          ijk.x += na;
00438       ijk.y = fmod(ijk.y,nb);
00439       if (ijk.y < 0)
00440          ijk.y += nb;
00441       ijk.z = fmod(ijk.z,nc);
00442       if (ijk.z < 0)
00443          ijk.z += nc;
00444       
00445       i = celli[(int)(floor(ijk.x))];
00446       j = cellj[(int)(floor(ijk.y))];
00447       k = cellk[(int)(floor(ijk.z))];
00448       return;
00449    }
00450    
00451    // get_xyz(const Vector3D ijk) const
00452    // Convert the grid coordinates ijk to real coordinates
00453    // The grid coordinates may have a decimal component, indicating
00454    // where it falls within the grid block.
00455    //
00456    Vector3D get_xyz(const Vector3D ijk) const
00457    {
00458       Vector3D v = grida * ijk.x / na + gridb * ijk.y / nb + gridc * ijk.z / nc 
00459       + origin;
00460       return v;
00461    }
00462    
00463    // get_xyz(const int i, const int j, const int k) const
00464    // Convert the grid coordinates ijk to real coordinates
00465    // The real coordinates returned are for the center of the grid block
00466    //
00467    Vector3D get_xyz(const int i, const int j, const int k) const 
00468    {
00469       const double id = (i+0.5) / na;
00470       const double jd = (j+0.5) / nb;
00471       const double kd = (k+0.5) / nc;
00472       
00473       Vector3D v = grida * id + gridb * jd + gridc * kd + origin;
00474       return v;
00475    }
00476    
00477    // get_cell(const int i, const int j, const int k) const
00478    // Convert the grid coordinates ijk to real coordinates
00479    // The real coordinates returned are for the center of the grid block
00480    //
00481    const GridList* get_cell(const int i, const int j, const int k) const
00482    {
00483       if (i>=0 && i<cna && j>=0 && j<cnb && k>=0 && k<cnc) {
00484          return cellgrid[((i*cnb)+j)*cnc+k]->get_list();
00485       } 
00486       return 0;
00487    }
00488    
00489    Vector3D wrap_xyz(const Vector3D xyz) const
00490    {
00491       Vector3D gridcoord = get_ijk(xyz);
00492       gridcoord.x = fmod(gridcoord.x,na);
00493       if (gridcoord.x < 0)
00494          gridcoord.x += na;
00495       gridcoord.y = fmod(gridcoord.y,nb);
00496       if (gridcoord.y < 0)
00497          gridcoord.y += nb;
00498       gridcoord.z = fmod(gridcoord.z,nc);
00499       if (gridcoord.z < 0)
00500          gridcoord.z += nc;
00501       
00502       const Vector3D ret = get_xyz(gridcoord);
00503       //    if (gc2.x != gridcoord.x || gc2.y != gridcoord.y || gc2.z != gridcoord.z) {
00504       //      printf("gc2 %f %f %f grid %f %f %f ret %f %f %f\n",
00505       //        gc2.x,gc2.y,gc2.z,gridcoord.x,
00506       //        gridcoord.y,gridcoord.z,ret.x,ret.y,ret.z);
00507       //    }
00508       return ret;
00509    }
00510    
00511    int get_na() const { return na; }
00512    int get_nb() const { return nb; }
00513    int get_nc() const { return nc; }
00514    
00515    int get_cna() const { return cna; }
00516    int get_cnb() const { return cnb; }
00517    int get_cnc() const { return cnc; }
00518    
00519    Vector3D findAtomBox(const double radius) const;
00520    void build(const MolData* mol_data);
00521    void store(char* fname);
00522    void print();
00523    static int* scramble_array(const int n);
00524    
00525 private:
00526       void Inverse3(const double m[3][3], double inv[3][3]) const;
00527    double Cofactor3(const double mat[3][3],const int i, const int j) const;
00528    
00529    
00530    Vector3D origin;
00531    Vector3D grida;
00532    Vector3D gridb;
00533    Vector3D gridc;
00534    double inv_abc[3][3];
00535    int na;
00536    int nb;
00537    int nc;
00538    
00539    int cna;
00540    int cnb;
00541    int cnc;
00542    AtomGridCell** cellgrid;
00543    int* celli;
00544    int* cellj;
00545    int* cellk;
00546    bool periodic;
00547 };
00548 
00549 
00550 MolData* MolData::readFile(const char* fname, double gridspacing, 
00551                            double radius, double threshold) {
00552   MolData *data = new MolData;
00553   data->gridspacing = gridspacing;
00554   data->radius = radius;
00555   data->threshold = threshold;
00556   data->actual_radius = radius+threshold;
00557    
00558    
00559   FILE* inpf = fopen(fname,"r");
00560   if (inpf == NULL) {
00561     printf("Couldn't open file %s\n",fname);
00562     exit(-1);
00563   }
00564    
00565   int n_read=0;
00566   n_read += fscanf(inpf,"%lf %lf %lf",
00567                    &(data->origin.x),&(data->origin.y),&(data->origin.z)); 
00568   n_read += fscanf(inpf,"%lf %lf %lf",
00569                    &(data->basisa.x),&(data->basisa.y),&(data->basisa.z)); 
00570   n_read += fscanf(inpf,"%lf %lf %lf",
00571                    &(data->basisb.x),&(data->basisb.y),&(data->basisb.z)); 
00572   n_read += fscanf(inpf,"%lf %lf %lf",
00573                    &(data->basisc.x),&(data->basisc.y),&(data->basisc.z));
00574   n_read += fscanf(inpf,"%d",&(data->count));
00575    
00576   if (n_read != 13) {
00577     printf("Error reading header (%d)",n_read);
00578     exit(-2);
00579   }
00580    
00581   data->coords = new AtomData[data->count];
00582   int i;
00583   for (i=0; i<data->count; i++) {
00584     n_read = fscanf(inpf,"%d %lf %lf %lf",
00585                     &(data->coords[i].index),
00586                     &(data->coords[i].coord.x),
00587                     &(data->coords[i].coord.y),
00588                     &(data->coords[i].coord.z));
00589     if (n_read != 4) {
00590        printf("Error reading atom %d (%d)",i,n_read);
00591        exit(-2);
00592     }
00593   }
00594   fclose(inpf);
00595   return data;
00596 }
00597 
00598 
00599 MolData* MolData::readData(const int *indx, const float *coords, 
00600                            const int count,
00601                            const Vector3D o,
00602                            const Vector3D a,
00603                            const Vector3D b,
00604                            const Vector3D c,
00605                            const double gridspacing,
00606                            const double radius,
00607                            const double threshold,
00608                            const bool periodic)
00609 {
00610   MolData *data = new MolData;
00611   data->gridspacing = gridspacing;
00612   data->radius = radius;
00613   data->threshold = threshold;
00614   data->actual_radius = radius+threshold;
00615   data->periodic = periodic;
00616    
00617   data->origin=o;
00618   data->basisa = a;
00619   data->basisb = b;
00620   data->basisc = c;
00621   data->count = count;
00622         
00623   data->coords = new AtomData[data->count];
00624   int i;
00625   int ind=0;
00626   for (i=0; i<data->count; i++) {
00627     data->coords[i].index = indx[i];
00628     data->coords[i].coord.x = coords[ind];
00629     data->coords[i].coord.y = coords[ind+1];
00630     data->coords[i].coord.z = coords[ind+2];
00631     ind += 3;
00632   }
00633   return data;
00634 }
00635 
00636 
00637 void MolData::wrapCoords(const AtomGrid* grid) {
00638   int i;
00639   for(i=0; i < count; i++) {
00640     coords[i].coord = grid->wrap_xyz(coords[i].coord);
00641   }
00642 }
00643 
00644 
00645 int GridList::add(const int i, const int j, const int k) {
00646   if (indx >= n)
00647     return -1;
00648    
00649   const int i2 = shuffle_vec[indx];
00650   items[i2].x = i;
00651   items[i2].y = j;
00652   items[i2].z = k;
00653   indx++;
00654   return indx-1;
00655 }
00656 
00657 
00658 int AtomGridCell::set(const int i, const int j, const int k) {
00659   if (i<mina || i>maxa || j<minb || j>maxb || k<minc || k>maxc ) {
00660     return -1;
00661   }
00662   if (getlistcalled) {
00663     return -2;
00664   }
00665    
00666   if (noneset) {
00667     cellmap = new bool[na*nb*nc];
00668     int m;
00669     for (m=na*nb*nc-1; m >= 0; m--)
00670       cellmap[m] = false;
00671     noneset = false;
00672     // msgInfo << "Allocated cell " << i << "," << j << "," << k << sendmsg;
00673   }
00674    
00675   const int ii = i - mina;
00676   const int jj = j - minb;
00677   const int kk = k - minc;
00678   const int indx = ((ii*nb)+jj)*nc+kk;
00679   if (!cellmap[indx]) {
00680     count++;
00681   } 
00682   cellmap[indx] = true;
00683  
00684   return 0;
00685 }
00686 
00687 
00688 bool AtomGridCell::get(const int i, const int j, const int k) const {
00689   //  std::cout << "AtomGridCell::get " << i << "," <<mina << "," << maxa 
00690   //    << std::endl;
00691   
00692   if (i<mina || i>maxa || j<minb || j>maxb || k<minc || k>maxc )
00693     return false;
00694 
00695   if (noneset)
00696     return false;
00697    
00698   const int ii = i - mina;
00699   const int jj = j - minb;
00700   const int kk = k - minc;
00701   return cellmap[((ii*nb)+jj)*nc+kk];
00702 }
00703 
00704 
00705 const GridList* AtomGridCell::get_list() {
00706   if (!getlistcalled) {
00707     build_list();
00708   }
00709   return gridlist;
00710 }
00711 
00712 
00713 void AtomGridCell::build_list() {
00714   // Build a list containing only vacuum cells
00715   getlistcalled = true;
00716   int sz = na*nb*nc;
00717   bool* vac = new bool[sz];
00718   int i;
00719   for (i=0; i < sz; i++)
00720     vac[i] = false;
00721    
00722   int tot=0;
00723   int grid_idx = 0;
00724   for(i=mina; i <= maxa; i++) {
00725     // msgInfo << i << sendmsg;
00726     int j;
00727     for (j=minb; j <= maxb; j++) {
00728       // msgInfo << j << "," << minc << "-" << maxc << ":";
00729       int k;
00730       for (k=minc; k <= maxc; grid_idx++, k++) {
00731         if (!get(i,j,k)) {
00732           // Optimization:
00733           // If all neighbors are also vacuum, then don't add this one
00734           // since its an interior cell
00735           bool all_vac = true;
00736           if (i > mina && i < maxa 
00737               && j > minb && j < maxb
00738               && k > minc && k < maxc) {
00739             int ii;
00740             for (ii=-1; ii <= 1; ii++) {
00741               int jj;
00742               for (jj=-1; jj <= 1; jj++) {
00743                 int kk;
00744                 for (kk=-1; kk <= 1; kk++) {
00745                   if (get(i+ii,j+jj,k+kk)) {
00746                     // msgInfo << "Not storing " << i << "," << j << "," << k 
00747                     //         << sendmsg;
00748                     all_vac = false;
00749                     break;
00750                   }
00751                 }
00752                 if (!all_vac)
00753                   break;
00754               }
00755               if (!all_vac)
00756                 break;
00757             }
00758           } else {
00759             all_vac = false;
00760           }
00761           if (!all_vac) {
00762             // msgInfo << "Storing " << i << "," << j << "," << k << sendmsg;
00763             vac[grid_idx] = true;
00764             tot++;
00765             // msgInfo << "1";
00766           } // else std::cout << "2";
00767         } // else std::cout << "0";
00768       }
00769       // msgInfo << sendmsg;
00770     }
00771   }
00772   gridlist = new GridList(ag,tot);
00773   grid_idx = 0;
00774   int gi;
00775   for (gi=mina; gi <= maxa; gi++) {
00776     int gj;
00777     for (gj=minb; gj <= maxb; gj++) {
00778       int gk;
00779       for (gk=minc; gk <= maxc; grid_idx++, gk++) {
00780         if (vac[grid_idx]) {
00781           // msgInfo << "Adding " << i << "," << j << "," << k << sendmsg;
00782           gridlist->add(gi,gj,gk);
00783         }
00784       }
00785     }
00786   }
00787 
00788   //  msgInfo << "Gridlist storing " << tot << " of " << sz 
00789   //          << " elements noneset=" << noneset << " count=" << count 
00790   //          << sendmsg;
00791   delete [] vac;
00792 }
00793 
00794 
00795 AtomGrid::AtomGrid(const MolData* mdata) {
00796   // We'll make a grid with rectangular cells of the specified size
00797   // with a corner at the origin. This means the grid may not exactly
00798   // fit with the periodic box
00799   origin = mdata->origin;
00800   grida = mdata->basisa;
00801   gridb = mdata->basisb;
00802   gridc = mdata->basisc;
00803   na = (int)ceil(mdata->basisa.length() / mdata->gridspacing);
00804   nb = (int)ceil(mdata->basisb.length() / mdata->gridspacing);
00805   nc = (int)ceil(mdata->basisc.length() / mdata->gridspacing);
00806   periodic = mdata->periodic;
00807    
00808   // We'll make a grid with rectangular cells of the specified size
00809   // with a corner at the origin. This means the grid may not exactly
00810   // fit with the periodic box
00811   printf("Grid %d %d %d\n",na,nb,nc);
00812   printf("Origin %f %f %f\n",origin.x,origin.y,origin.z);
00813    
00814   double a[3][3];
00815    
00816   a[0][0] = grida.x;
00817   a[0][1] = gridb.x;
00818   a[0][2] = gridc.x;
00819   a[1][0] = grida.y;
00820   a[1][1] = gridb.y;
00821   a[1][2] = gridc.y;
00822   a[2][0] = grida.z;
00823   a[2][1] = gridb.z;
00824   a[2][2] = gridc.z;
00825   Inverse3(a,this->inv_abc);
00826    
00827   printf("a...\n");
00828   printf("%f %f %f\n",a[0][0],a[0][1],a[0][2]);
00829   printf("%f %f %f\n",a[1][0],a[1][1],a[1][2]);
00830   printf("%f %f %f\n",a[2][0],a[2][1],a[2][2]);
00831   
00832   printf("inv a...\n");
00833   printf("%f %f %f\n",inv_abc[0][0],inv_abc[0][1],inv_abc[0][2]);
00834   printf("%f %f %f\n",inv_abc[1][0],inv_abc[1][1],inv_abc[1][2]);
00835   printf("%f %f %f\n",inv_abc[2][0],inv_abc[2][1],inv_abc[2][2]);
00836    
00837   // Find how many cells to divide this into
00838   Vector3D cellsz = findAtomBox(mdata->getActualRadius());
00839   printf("cellsz=%f, %f, %f\n",cellsz.x,cellsz.y,cellsz.z);
00840    
00841   cna = (int)(floor(na/cellsz.x));
00842   cnb = (int)(floor(nb/cellsz.y));
00843   cnc = (int)(floor(nc/cellsz.z));
00844   if (cna==0) {
00845     cna = 1;
00846   }
00847   if (cnb==0) {
00848     cnb = 1;
00849   }
00850   if (cnc==0) {
00851     cnc = 1;
00852   }
00853    
00854   printf("na=%d, %d, %d\n",na,nb,nc);
00855   printf("cna=%d, %d, %d\n",cna,cnb,cnc);
00856    
00857   // Build the cell grid, dividing grid points roughly-evenly among the cells
00858   cellgrid = new AtomGridCell*[cna*cnb*cnc];
00859   // Storing the cell boundaries will let us find the cells quickly.
00860   celli = new int[na+1];
00861   celli[na] = cna;
00862   cellj = new int[nb+1];
00863   cellj[nb] = cnb;
00864   cellk = new int[nc+1];
00865   cellk[nc] = cnc;
00866    
00867   int* gridi = new int[cna+1];
00868   gridi[cna] = na;
00869   int* gridj = new int[cnb+1];
00870   gridj[cnb] = nb;
00871   int* gridk = new int[cnc+1];
00872   gridk[cnc] = nc;
00873    
00874   // Build the mapping from grid point to cells
00875   int mina = 0;
00876   int maxa;
00877   int i;
00878   for (i=0; i < cna; i++) {
00879     maxa = ((i+1) * na) / cna - 1;
00880     // msgInfo << "minmax=" << mina << "," << maxa << sendmsg;
00881     int ii;
00882     for (ii=mina; ii<= maxa; ii++) {
00883       celli[ii] = i;
00884       // msgInfo << "cell[" << ii << "]=" << celli[ii] << sendmsg;
00885     }
00886     gridi[i] = mina;
00887     int minb = 0;
00888     int maxb;
00889     int j;
00890     for(j=0; j < cnb; j++) {
00891       maxb = ((j+1) * nb) / cnb - 1;
00892       if (i==0) {
00893         for (int jj=minb; jj<= maxb; jj++)
00894           cellj[jj] = j;
00895       }
00896       gridj[j] = minb;
00897       int minc = 0;
00898       int maxc;
00899       int k;
00900       for (k=0; k < cnc; k++) {
00901         maxc = ((k+1) * nc) / cnc - 1;
00902         if (i==0 && j==0) {
00903           for (int kk=minc; kk<= maxc; kk++)
00904             cellk[kk] = k;
00905         }
00906         gridk[k] = minc;
00907         minc = maxc+1;
00908       }
00909       minb = maxb+1;
00910     }
00911     mina = maxa+1;
00912   }
00913    
00914   // Now create the cells
00915   int ci;
00916   for (ci=0; ci < cna; ci++) {
00917     int cj;
00918     for (cj=0; cj < cnb; cj++) {
00919       int ck;
00920       for (ck=0; ck < cnc; ck++) {
00921         // msgInfo << "cellgrid[" << i <<"," << j << "," << k << "," <<   ((i*cnb)+j)*cnc + k << " addr " << cellgrid << sendmsg;
00922         cellgrid[((ci*cnb)+cj) * cnc + ck] = 
00923           new AtomGridCell(this,gridi[ci],gridi[ci+1]-1,
00924                                 gridj[cj],gridj[cj+1]-1,
00925                                 gridk[ck],gridk[ck+1]-1);
00926       }
00927     }
00928   }
00929 
00930   delete[] gridi;
00931   delete[] gridj;
00932   delete[] gridk;
00933 }
00934 
00935 
00936 int* GridList::scramble_array(const int n) {
00937    return AtomGrid::scramble_array(n); 
00938 }
00939 
00940 
00941 int* AtomGrid::scramble_array(const int n) {
00942   if (n < 1) return 0;
00943    
00944   int* a = new int[n];
00945   a[0]=n-1;
00946   if (n == 1) return a;
00947   a[1] = 0;
00948   if (n == 2) return a;
00949    
00950   bool* b = new bool[n];
00951   int i;
00952   for (i=0; i<n; i++)
00953     b[i] = true;
00954    
00955   int howmany = 1;
00956   int incr = (n >> 1);
00957   // msgInfo << "n=" << n << " n>>=" << (n >> 1) << sendmsg;
00958    
00959   int nxt=2;
00960   while (nxt < n) {
00961     int i;
00962     for (i=0; i < howmany; i++) {
00963       //      std::cout << "n=" << n << " nxt=" << nxt << " i=" << i << " incr=" << incr << " howmany=" << howmany;
00964       int nxtval = a[i+1] + incr;
00965       //      std::cout << "  nxtval=" << nxtval;
00966       if (b[nxtval])  {
00967         a[nxt++] = nxtval;
00968         //        std::cout << " a[" << nxt-1 << "]=" << a[nxt-1];
00969         b[nxtval] = false;
00970       } 
00971       //      std::cout << std::endl;
00972       if (nxt > n-1)
00973         break;
00974     }
00975     howmany <<= 1;
00976     incr >>= 1;
00977     if (incr == 0)
00978       break;
00979   }
00980 
00981   for (i=1; i < n; i++) {
00982     if (nxt >= n) break;
00983     if (b[i])
00984       a[nxt++] = i;
00985   }
00986 
00987   delete [] b;
00988   return a;
00989 }
00990 
00991 
00992 Vector3D AtomGrid::findAtomBox(const double radius) const {
00993   // Find the bounding box for the atom radius
00994   Vector3D cornerx[8];
00995   Vector3D dx, dy, dz;
00996    
00997   const Vector3D o = origin;
00998   dx.y = dx.z = dy.x = dy.z = dz.x = dz.y = 0;
00999   dx.x = dy.y = dz.z = radius;
01000   cornerx[0] = o;
01001   cornerx[1] = o + dx;
01002   cornerx[2] = o + dy;
01003   cornerx[3] = o + dz;
01004   cornerx[4] = cornerx[1] + dy;
01005   cornerx[5] = cornerx[1] + dz;
01006   cornerx[6] = cornerx[2] + dz;
01007   cornerx[7] = cornerx[4] + dz;  
01008    
01009   // Find the max i,j,k for those corners
01010   Vector3D ijk = get_ijk(cornerx[0]);
01011   double dimax = ijk.x;
01012   double djmax = ijk.y;
01013   double dkmax = ijk.z;
01014   
01015   int i;
01016   for (i=1; i<8; i++) {
01017     ijk = get_ijk(cornerx[i]);
01018     if (ijk.x > dimax)
01019       dimax = ijk.x;
01020     if (ijk.y > djmax)
01021       djmax = ijk.y;
01022     if (ijk.z > dkmax)
01023       dkmax = ijk.z;
01024   }
01025   Vector3D vec;
01026    
01027   vec.x = dimax;
01028   vec.y = djmax;
01029   vec.z = dkmax;
01030    
01031   return vec;
01032 }
01033 
01034 
01035 void AtomGrid::build(const MolData* mdata) {
01036   // Fill the grid
01037   // Find the bounding box for the atom radius
01038   Vector3D bb = findAtomBox(mdata->getRadius());
01039    
01040   const int imax = (int)ceil(bb.x);
01041   const int jmax = (int)ceil(bb.y);
01042   const int kmax = (int)ceil(bb.z);
01043    
01044   printf("Box=%d, %d, %d\n",imax,jmax,kmax);
01045   const double rsq = mdata->getRadius() * mdata->getRadius();
01046   int atomi;
01047   for (atomi=0; atomi<mdata->count; atomi++) {
01048     Vector3D atom = mdata->coords[atomi].coord;
01049 //    msgInfo << "Checking atom " << atomi 
01050 //             << "(" << atom.x << "," << atom.y << "," << atom.z << sendmsg;
01051     Vector3D atomijk = get_ijk(atom);
01052     const int ai = (int)atomijk.x;
01053     const int aj = (int)atomijk.y;
01054     const int ak = (int)atomijk.z;
01055     int i;
01056     for (i=ai-imax; i <= ai+imax; i++) {
01057       int ii = i;
01058       if (i < 0) {
01059         if (periodic) {
01060           ii += na;
01061         } else continue;
01062       } else if (i >= na) {
01063         if (periodic) {
01064           ii -= na;
01065         } else continue;
01066       }
01067 
01068       int j;
01069       for(j=aj-jmax; j <= aj+jmax; j++) {
01070         int jj = j;
01071         if (j < 0) {
01072           if (periodic) {
01073             jj += nb;
01074           } else continue;
01075         } else if (j >= nb) {
01076           if (periodic) {
01077             jj -= nb;
01078           } else continue;
01079         }
01080   
01081         // Do the positive half of the loop. Stop once we get to
01082         // the first cell that is outside the boundary
01083         int k;
01084         for (k=0; k <= kmax; k++) {
01085           int k_no_wrap = ak+k;
01086           int kk = k_no_wrap;
01087           if (k_no_wrap < 0) {
01088             if (periodic) {
01089               kk += nc;
01090             } else continue;
01091           } else if (k_no_wrap >= nc) {
01092             if (periodic) {
01093               kk -= nc;
01094             } else continue;
01095           }
01096            
01097           // If cell is already filled, just continue
01098           if (get(ii,jj,kk))
01099             continue;
01100                
01101           const Vector3D v = get_xyz(i,j,k_no_wrap);
01102           const Vector3D dv = v - atom;
01103           if (dv.length2() <= rsq) {
01104             set(ii,jj,kk);
01105           } else {
01106             break;
01107           }
01108         }
01109 
01110         for (k=1; k <= kmax; k++) {
01111           int k_no_wrap = ak-k;
01112           int kk = k_no_wrap;
01113           if (k_no_wrap < 0) {
01114             if (periodic) {
01115               kk += nc;
01116             } else continue;
01117           } else if (k_no_wrap >= nc) {
01118             if (periodic) {
01119               kk -= nc;
01120             } else continue;
01121           }
01122           
01123           // If cell is already filled, just continue
01124           if (get(ii,jj,kk))
01125             continue;
01126           
01127           const Vector3D v = get_xyz(i,j,k_no_wrap);
01128           const Vector3D dv = v - atom;
01129           if (dv.length2() <= rsq) {
01130             set(ii,jj,kk);
01131           } else {
01132             break;
01133           }
01134         }
01135       }
01136     }
01137   }
01138   return;
01139 }
01140 
01141 
01142 void AtomGrid::print() {
01143   int k;
01144   for (k=0; k<nc; k++) {
01145     int i;
01146     for (i=0; i<na; i++) {
01147       printf("%d:",i);
01148       int j;
01149       for (j=0; j<nb; j++) {
01150         if (get(i,j,k))
01151           printf("1");
01152         else printf("0");
01153       }
01154       printf("\n");
01155     }
01156     printf("%d:-------------------------------------------------------\n",k);
01157   }
01158 }
01159 
01160 
01161 void AtomGrid::store(char* fname) {
01162   FILE* outf = fopen(fname,"w");
01163   if (outf == NULL) {
01164     printf("Couldn't open file %s\n",fname);
01165     exit(-1);
01166   }
01167    
01168   int k;
01169   for (k=0; k<nc; k++) {
01170     int i;
01171     for (i=0; i<na; i++) {
01172       int j;
01173       for (j=0; j<nb; j++) {
01174         if (!get(i,j,k)) {
01175           Vector3D v = get_xyz(i,j,k);
01176           if (fprintf(outf,"%f %f %f\n",v.x,v.y,v.z) < 0) {
01177             printf("Store grid error %s\n",fname);
01178             exit(-2);
01179           }
01180         }
01181       }
01182     }
01183   }
01184    
01185   fclose(outf);
01186   return;
01187 }
01188 
01189 
01190 void AtomGrid::Inverse3(const double m[3][3], double inv[3][3]) const {
01191   // Get adjoint matrix : transpose of cofactor matrix
01192   int i,j;
01193   for (i=0; i<3; i++)
01194     for (j=0; j<3; j++)
01195       inv[i][j] = Cofactor3(m,i,j);
01196   
01197   // Now divide by the determinant
01198   const double det = m[0][0] * inv[0][0] 
01199                    + m[0][1] * inv[1][0]
01200                    + m[0][2] * inv[2][0];
01201    
01202   //    if (det == 0) {
01203   //      error "Non-invertible"
01204   //    }
01205    
01206   for (i=0; i<3; i++)
01207     for (j=0; j<3; j++)
01208       inv[i][j] /= det;
01209   
01210   return;
01211 }
01212 
01213 
01214 double AtomGrid::Cofactor3(const double mat[3][3],
01215                            const int i, const int j) const {
01216   int cols[3][2] = { {1,2}, {0,2}, {0,1}};
01217    
01218   const int row1 = cols[j][0];
01219   const int row2 = cols[j][1];
01220   
01221   const int col1 = cols[i][0];
01222   const int col2 = cols[i][1];
01223   
01224   const double a = mat[row1][col1];
01225   const double b = mat[row1][col2];
01226   const double c = mat[row2][col1];
01227   const double d = mat[row2][col2];
01228   
01229   const double det = a*d-b*c;
01230   if (((i+j) % 2) != 0)
01231     return -det;
01232   else
01233     return det;
01234 }
01235 
01236 
01237 AtomList* MolData::findAtomsGrid(const AtomGrid* grid) {
01238   // Scan atoms and check adjacent cells
01239   AtomList *al = new AtomList(count);
01240   const double r = getActualRadius();
01241   const double rsq = r*r;
01242   const Vector3D atom_box = grid->findAtomBox(radius);
01243 
01244   int atomi;
01245   for (atomi=0; atomi<count; atomi++) {
01246     Vector3D atom = coords[atomi].coord;
01247 //    msgInfo << "Checking atom " << atomi 
01248 //            << "(" << atom.x << "," << atom.y << "," << atom.z << sendmsg;
01249     int i0, j0, k0;
01250     bool found = false;
01251       
01252     grid->get_cijk(atom,i0,j0,k0);
01253     bool cell_outside = false;
01254 
01255     int di;
01256     for (di=-1; !found && di <= 1; di++) {
01257       int ci = i0 + di;
01258       int cna = grid->get_cna();
01259       int wrapi = 0;
01260       if (ci<0) {
01261         ci += cna;
01262         wrapi = grid->get_na();
01263         if (!periodic)
01264           cell_outside = true;
01265       } else if (ci >= cna) {
01266         ci -= cna;
01267         wrapi = -grid->get_na();
01268         if (!periodic)
01269           cell_outside = true;
01270       }
01271 
01272       int dj;
01273       for (dj=-1; !found && dj <= 1; dj++) {
01274         int cj = j0 + dj;
01275         int cnb = grid->get_cnb();
01276         int wrapj = 0;
01277         if (cj<0) {
01278           cj += cnb;
01279           wrapj = grid->get_nb();
01280           if (!periodic)
01281             cell_outside = true;
01282         } else if (cj >= cnb) {
01283           cj -= cnb;
01284           wrapj = grid->get_nb();
01285           if (!periodic)
01286             cell_outside = true;
01287         }
01288 
01289         int dk;
01290         for (dk=-1; !found && dk <= 1; dk++) {
01291           int ck = k0 + dk;
01292           int cnc = grid->get_cnc();
01293           int wrapk = 0;
01294           if (ck<0) {
01295             ck += cnc;
01296             wrapk = grid->get_nc();
01297             if (!periodic)
01298               cell_outside = true;
01299           } else if (ck >= cnc) {
01300             ck -= cnc;
01301             wrapk = -grid->get_nc();
01302             if (!periodic)
01303               cell_outside = true;
01304           }
01305 
01306           // If the cell is outside the box, then there must be a
01307           // vacuum there
01308           if (!periodic && cell_outside) {
01309             const Vector3D atom_grid = grid->get_ijk(atom);
01310             const Vector3D corner0 = atom_grid - atom_box;
01311             bool addit = false;
01312             if (corner0.x < 0 || corner0.y < 0 || corner0.z < 0) {
01313               addit = true;
01314             } else {
01315               const Vector3D corner1 = atom_grid + atom_box;
01316               if (corner1.x >= grid->get_na() 
01317                   || corner1.y >= grid->get_nb() 
01318                   || corner1.z >= grid->get_nc() )
01319                 addit = true;
01320             }
01321             if (addit) {
01322               al->add(coords[atomi].index);
01323               found = true;
01324             }
01325             cell_outside = false;
01326           }
01327 
01328           const GridList* list = grid->get_cell(ci,cj,ck);
01329 //        msgInfo << "Found " << list->getNum() 
01330 //                << " atoms in cell " 
01331 //                << ci << "," << cj << "," << ck << sendmsg;
01332           int gli;
01333           for (gli=0; !found && gli < list->getNum(); gli++) {
01334 //          XXX: this does not handle get() receiving an illegal argument.
01335             Vector3D gridpt = list->get(gli);
01336             gridpt.x += wrapi;
01337             gridpt.y += wrapj;
01338             gridpt.z += wrapk;
01339             const Vector3D v = grid->get_xyz(gridpt);
01340             const Vector3D dv = v - atom;
01341             if (dv.length2() <= rsq) {
01342               al->add(coords[atomi].index);
01343               found = true;
01344 //            msgInfo << "Found vac " << gli << sendmsg;
01345             }
01346           }
01347         }
01348       }
01349     }
01350   }
01351 
01352   return al;
01353 }
01354 
01355 
01356 int msmain(int argc,char *argv[]) {
01357   if (argc != 6) {
01358     printf("Usage: %s grid-res vac-radius selection-dist coord-file index-file\n", argv[0]);
01359      return 0;
01360   }
01361    
01362   double res = strtod(argv[1],NULL); // Max grid resolution
01363   double radius = strtod(argv[2],NULL); // vacuum-sphere/atom mask radius
01364   double select_dist = strtod(argv[3],NULL); // selection distance
01365    
01366   printf("Reading\n");
01367   MolData* mol_data = MolData::readFile(argv[4],res,radius,select_dist);
01368   printf("Building grid\n");
01369   AtomGrid* atom_grid = new AtomGrid(mol_data);
01370   mol_data->wrapCoords(atom_grid);
01371   atom_grid->build(mol_data);
01372    
01373   //  atom_grid->print();
01374   //  atom_grid->store("grid.dat");
01375    
01376   printf("Searching grid\n");
01377   AtomList* atom_list = mol_data->findAtomsGrid(atom_grid);
01378   printf("Storing\n");
01379   atom_list->storeFile(argv[5]);
01380    
01381   delete atom_list;
01382   delete atom_grid;
01383   delete mol_data;
01384   
01385   return 0;
01386 }
01387 
01388 
01389 void measure_surface_int(const double res, const double radius, 
01390                          const double select_dist,
01391                          const bool periodic,
01392                          const double a, const double b, const double c, 
01393                          const double alpha, const double beta, const double gamma,
01394                          int *indices, float *coords, int count,
01395                          int **surf, int *n_surf
01396                          ) 
01397 {
01398   if (periodic)
01399     printf("System periodic\n");
01400   else
01401     printf("System not periodic\n");
01402 
01403   printf("Reading\n");
01404   Vector3D uo;
01405   Vector3D ua;
01406   Vector3D ub;
01407   Vector3D uc;
01408    
01409   // Set origin to center of box
01410   double xmin = coords[0];
01411   double xmax = coords[0];
01412   double ymin = coords[1];
01413   double ymax = coords[1];
01414   double zmin = coords[2];
01415   double zmax = coords[2];
01416   int ind;
01417   int i;
01418   for(i=1, ind=3; i < count; i++, ind += 3) {
01419     if (coords[ind] < xmin)
01420       xmin = coords[ind];
01421     if (coords[ind] > xmax)
01422       xmax = coords[ind];
01423     if (coords[ind+1] < ymin)
01424       ymin = coords[ind+1];
01425     if (coords[ind+1] > ymax)
01426       ymax = coords[ind+1];
01427     if (coords[ind+2] < zmin)
01428       zmin = coords[ind+2];
01429     if (coords[ind+2] > zmax)
01430       zmax = coords[ind+2];
01431   }
01432   uo.x = 0.5*(xmin+xmax);
01433   uo.y = 0.5*(ymin+ymax);
01434   uo.z = 0.5*(zmin+zmax);
01435    
01436   double cosbc = cos(DEGTORAD(alpha));
01437   double cosac = cos(DEGTORAD(beta));
01438   double sinab, cosab;
01439   sincos(DEGTORAD(gamma), &sinab, &cosab);
01440 
01441   ua.x = a; ua.y = 0; ua.z = 0;
01442   ub.x = b * cosab;
01443   ub.y = b * sinab;
01444   ub.z = 0;
01445 
01446   // If sinAB is zero, then we can't determine C uniquely since it's defined
01447   // in terms of the angle between A and B.
01448   if (sinab > 0) {
01449     uc.x = cosac;
01450     uc.y = (cosbc - cosac * cosab) / sinab;
01451     uc.z = sqrt(1.0 - uc.x*uc.x - uc.y*uc.y);
01452   }
01453   uc = c * uc;
01454   // Now move origin to the bottom corner
01455   uo = uo - (0.5 * (ua + ub + uc));
01456    
01457   MolData* mol_data = MolData::readData(indices, coords, count,
01458                                         uo, ua, ub, uc,
01459                                         res,radius,select_dist,periodic);
01460   printf("Building grid\n");
01461   AtomGrid* atom_grid = new AtomGrid(mol_data);
01462   mol_data->wrapCoords(atom_grid);
01463   atom_grid->build(mol_data);
01464    
01465   //   atom_grid->print();
01466   //   atom_grid->store("grid.dat");
01467    
01468   printf("Searching grid\n");
01469   AtomList* atom_list = mol_data->findAtomsGrid(atom_grid);
01470   printf("Storing\n");
01471   atom_list->storeData(surf,n_surf);
01472    
01473   delete atom_list;
01474   delete atom_grid;
01475   delete mol_data;
01476      
01477   return;
01478 }
01479 
01480 
01481 int measure_surface(AtomSel *sel, MoleculeList *mlist,
01482                    const float *framepos,
01483                    const double gridsz, 
01484                    const double radius,
01485                    const double depth,
01486                    int **surface, int *n_surf ) {
01487   int i;
01488   
01489   if (!sel)                     return MEASURE_ERR_NOSEL;
01490   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
01491   
01492   // compute the axis aligned aligned bounding box for the selected atoms
01493   int *sel_ind = new int[sel->selected];
01494   float *coords = new float[sel->selected*3];
01495   int j=0;
01496   int k=0;
01497   int ind=sel->firstsel*3;
01498   for(i=sel->firstsel; i<=sel->lastsel; i++) {
01499     if (sel->on[i]) {
01500       sel_ind[j++] = i;
01501       coords[k++] = framepos[ind];
01502       coords[k++] = framepos[ind+1];
01503       coords[k++] = framepos[ind+2];
01504     }
01505     ind += 3;
01506   }
01507   
01508   Timestep *ts = sel->timestep(mlist);
01509   double a = ts->a_length;
01510   double b = ts->b_length;
01511   double c = ts->c_length;
01512   double alpha = ts->alpha;
01513   double beta = ts->beta;
01514   double gamma = ts->gamma;
01515   bool periodic = true;
01516   if (a == 0 || b == 0 || c == 0 || alpha == 0 || beta == 0 || gamma == 0) {
01517     periodic = false;
01518     float min_coord[3];
01519     float max_coord[3];
01520     measure_minmax(sel->num_atoms, sel->on, framepos, NULL, min_coord, max_coord);
01521     a = max_coord[0] - min_coord[0];
01522     b = max_coord[1] - min_coord[1];
01523     c = max_coord[2] - min_coord[2];
01524     alpha = beta = gamma = 90;
01525   }
01526         
01527   printf("Periodic: %f %f %f %f %f %f\n",a,b,c,alpha,beta,gamma);
01528   
01529   measure_surface_int(gridsz, radius, depth, periodic,
01530                       a,b,c,alpha,beta,gamma,
01531                       sel_ind,coords,sel->selected,
01532                       surface, n_surf);
01533   
01534   printf("Done!\n");
01535   
01536   return MEASURE_NOERR;
01537 }

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