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

MeasurePBC.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: MeasurePBC.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.18 $       $Date: 2019/01/17 21:21:00 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code to measure atom distances, angles, dihedrals, etc, 
00019  *   accounting for periodic boundary conditions
00020  ***************************************************************************/
00021 
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include "Measure.h"
00026 #include "AtomSel.h"
00027 #include "Matrix4.h"
00028 #include "utilities.h"
00029 #include "SpatialSearch.h"  // for find_within()
00030 #include "MoleculeList.h"
00031 #include "Inform.h"
00032 #include "Timestep.h"
00033 #include "VMDApp.h"
00034 
00035 //
00036 // Find an orthogonal basis R^3 with ob1=b1     
00037 //
00038 
00039 void orthonormal_basis(const float b[9], float e[9]) {
00040   float ob[3*3];
00041   vec_copy(ob+0, b+0);
00042   vec_copy(e+0, ob+0);
00043   vec_normalize(e+0);
00044   vec_triad(ob+3, b+3, -dot_prod(e+0, b+3), e+0);
00045   vec_copy(e+3, ob+3);
00046   vec_normalize(e+3);
00047   vec_triad(ob+6,  b+6, -dot_prod(e+0, b+6), e+0);
00048   vec_triad(ob+6, ob+6, -dot_prod(e+3, b+6), e+3);
00049   vec_copy(e+6, ob+6);
00050   vec_normalize(e+6);
00051 }
00052 
00053 //
00054 // Returns basis vectors in coordinates of an orthonormal 
00055 // basis obase.                                        
00056 //
00057 
00058 void basis_change(const float *base, const float *obase, float *newcoor, int n) {
00059   int i, j;
00060   for (i=0; i<n; i++) {
00061     for (j=0; j<n; j++) {
00062       newcoor[n*i+j] = dot_prod(&base[n*j], &obase[n*i]);
00063     }
00064   }
00065 }
00066 
00067 // Compute matrix that transforms coordinates from an arbitrary PBC cell 
00068 // into an orthonormal unitcell. Since the cell origin is not stored by VMD
00069 // you have to specify it.
00070 int measure_pbc2onc(MoleculeList *mlist, int molid, int frame, const float origin[3], Matrix4 &transform) {
00071   int orig_ts, max_ts;
00072 
00073   Molecule *mol = mlist->mol_from_id(molid);
00074   if( !mol )
00075     return MEASURE_ERR_NOMOLECULE;
00076  
00077   // get current frame number and make sure there are frames
00078   if((orig_ts = mol->frame()) < 0)
00079     return MEASURE_ERR_NOFRAMES;
00080   
00081   // get the max frame number and determine frame range
00082   max_ts = mol->numframes()-1;
00083   if (frame==-2)  frame = orig_ts;
00084   else if (frame>max_ts || frame==-1) frame = max_ts;
00085 
00086   Timestep *ts = mol->get_frame(frame);
00087 
00088   Matrix4 AA, BB, CC;
00089   ts->get_transforms(AA, BB, CC);
00090 
00091   // Construct the cell spanning vectors
00092   float cell[9];
00093   cell[0] = AA.mat[12];
00094   cell[1] = AA.mat[13];
00095   cell[2] = AA.mat[14];
00096   cell[3] = BB.mat[12];
00097   cell[4] = BB.mat[13];
00098   cell[5] = BB.mat[14];
00099   cell[6] = CC.mat[12];
00100   cell[7] = CC.mat[13];
00101   cell[8] = CC.mat[14];
00102 
00103   get_transform_to_orthonormal_cell(cell, origin, transform);
00104 
00105   return MEASURE_NOERR;
00106 }
00107 
00108 // Compute matrix that transforms coordinates from an arbitrary cell 
00109 // into an orthonormal unitcell. Since the origin is not stored by VMD
00110 // you have to specify it.
00111 // This is the lowlevel backend of measure_pbc2onc().
00112 
00113 // Here is a 2D example:
00114 // A and B are the are the displacement vectors which are needed to create 
00115 // the neighboring images. The parallelogram denotes the PBC cell with the origin O at its center.
00116 // The sqare to the right indicates the orthonormal unit cell i.e. the area into which the atoms 
00117 // will be wrapped by transformation T.
00118 //
00119 //                  + B                                        
00120 //                 /                              ^ B'         
00121 //       _________/________                       |            
00122 //      /        /        /                   +---|---+        
00123 //     /        /        /              T     |   |   |        
00124 //    /        O--------/-------> A   ====>   |   O---|--> A'  
00125 //   /                 /                      |       |        
00126 //  /_________________/                       +-------+        
00127   
00128 //  A  = displacement vector along X-axis with length a
00129 //  B  = displacement vector in XY-plane with length b
00130 //  A' = displacement vector along X-axis with length 1
00131 //  B' = displacement vector along Y-axis with length 1
00132 //  O = origin of the unit cell 
00133 
00134 void get_transform_to_orthonormal_cell(const float *cell, const float *center, Matrix4 &transform) {
00135   // Orthogonalize system:
00136   // Find an orthonormal basis of the cell (in cartesian coords).
00137   // If the cell vectors from VMD/NAMD are used this should actually always
00138   // return the identity matrix due to the way the cell vectors A, B and C
00139   // are defined (i.e. A || x; B lies in the x,y-plane; A, B, C form a right
00140   // hand system).
00141   float obase[3*3];
00142   orthonormal_basis(cell, obase);
00143 
00144   // Get orthonormal base in cartesian coordinates (it is the inverse of the
00145   // obase->cartesian transformation):
00146   float identity[3*3] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
00147   float obase_cartcoor[3*3];
00148   basis_change(obase, identity, obase_cartcoor, 3);
00149 
00150 
00151   // Transform 3x3 into 4x4 matrix:
00152   Matrix4 obase2cartinv;
00153   trans_from_rotate(obase_cartcoor, &obase2cartinv);
00154 
00155   // This is the matrix for the obase->cartesian transformation:
00156   Matrix4 obase2cart = obase2cartinv;
00157   obase2cart.inverse();
00158 
00159   // Get coordinates of cell in terms of obase
00160   float m[3*3]; 
00161   basis_change(cell, obase, m, 3);
00162   Matrix4 rotmat;
00163   trans_from_rotate(m, &rotmat);
00164   rotmat.inverse();
00165 
00166   
00167   // Actually we have:
00168   // transform = translation * obase2cart * obase2cartinv * rotmat * obase2cart
00169   //                           `------------v------------'
00170   //                                       = 1
00171   transform = obase2cart;
00172   transform.multmatrix(rotmat); // pre-multiplication
00173 
00174   // Finally we need to apply the translation of the origin
00175   float origin[3];
00176   vec_copy(origin, center);
00177   vec_scaled_add(origin, -0.5, &cell[0]);
00178   vec_scaled_add(origin, -0.5, &cell[3]);
00179   vec_scaled_add(origin, -0.5, &cell[6]);
00180   vec_negate(origin, origin);
00181   //printf("origin={%g %g %g}\n", origin[0], origin[1], origin[2]);
00182   transform.translate(origin);
00183 }
00184 
00185 
00186 // Get the array of coordinates of pbc image atoms for the specified selection.
00187 // The cutoff vector defines the region surrounding the pbc cell for which image 
00188 // atoms shall be constructed ({6 8 0} means 6 Angstrom for the direction of A,
00189 // 8 for B and no images in the C direction). The indexmap_array relates each
00190 // image atom to its corresponding main atom.
00191 // In case the molecule was aligned you can supply the alignment matrix which
00192 // is then used to correct for the rotation and shift of the pbc cell.
00193 // Since the pbc cell center is not stored in Timestep it must be provided.
00194 
00195 // The algorithm transforms the unitcell so that the unitcell minus the cutoff fits
00196 // into an orthonormal cell. Now the atoms in the rim can be easily identified and
00197 // wrapped into the neigboring cells. This works only if the largest cutoff
00198 // dimension is smaller than half of the smallest cell dimension. Otherwise a
00199 // slower algorithm is used that wraps each atom into all 26 neighboring cells
00200 // and checks if the image lies within cutoff.
00201 //
00202 //       ________________          
00203 //      / ____________  /          +---------+
00204 //     / /           / /           | +-----+ |
00205 //    / /   core    / /   ---->    | |     |_|___orthonormal_cell
00206 //   / /___________/ /    <----    | |     | |
00207 //  /_______________/              | +-----+ |___rim
00208 //                                 +---------+
00209 //
00210 // Alternatively one can specify a rectangular bounding box into which atoms
00211 // shall be wrapped. It is specified in form of minmax coordinates through
00212 // parameter *box. I.e. coordinates are produced for pbc image atom that lie
00213 // inside the box but outside the central unit cell. This feature can be used
00214 // for instance to retrieve coordinates of the minmax box of a selection when
00215 // the box boundaries exceed the unit cell.
00216 //
00217 // If a selection is provided (sel!=NULL) we return only coordinates that are
00218 // within the cutoff distance of that selection:
00219 
00220 // The results are provided in form of an array of 'extended' coordinates,
00221 // i.e. the coordinates of the requested region that don't lie in the central
00222 // unit cell. In order to identify these coordinates with the respective atoms
00223 // in the central cell an index map is also provided.
00224 
00225 // TODO:
00226 // * Allow requesting specific neighbor cells.
00227 
00228 int measure_pbc_neighbors(MoleculeList *mlist, AtomSel *sel, int molid,
00229                           int frame, const Matrix4 *alignment,
00230                           const float *center, const float *cutoff, const float *box,
00231                           ResizeArray<float> *extcoord_array,
00232                           ResizeArray<int> *indexmap_array) {
00233   int orig_ts, max_ts;
00234   if (!box && !cutoff[0] && !cutoff[1] && !cutoff[2]) return MEASURE_NOERR;
00235 
00236   Molecule *mol = mlist->mol_from_id(molid);
00237   if( !mol )
00238     return MEASURE_ERR_NOMOLECULE;
00239  
00240   // get current frame number and make sure there are frames
00241   if((orig_ts = mol->frame()) < 0)
00242     return MEASURE_ERR_NOFRAMES;
00243   
00244   // get the max frame number and determine current frame
00245   max_ts = mol->numframes()-1;
00246   if (frame==-2)  frame = orig_ts;
00247   else if (frame>max_ts || frame==-1) frame = max_ts;
00248 
00249   Timestep *ts = mol->get_frame(frame);
00250   if (!ts) return MEASURE_ERR_NOMOLECULE;
00251 
00252   // Get the displacement vectors (in form of translation matrices)
00253   Matrix4 Tpbc[3][2];
00254   ts->get_transforms(Tpbc[0][1], Tpbc[1][1], Tpbc[2][1]);
00255 
00256   // Assign the negative cell translation vectors
00257   Tpbc[0][0] = Tpbc[0][1];
00258   Tpbc[1][0] = Tpbc[1][1];
00259   Tpbc[2][0] = Tpbc[2][1];
00260   Tpbc[0][0].inverse();
00261   Tpbc[1][0].inverse();
00262   Tpbc[2][0].inverse();
00263 
00264   // Construct the cell spanning vectors
00265   float cell[9];
00266   cell[0] = Tpbc[0][1].mat[12];
00267   cell[1] = Tpbc[0][1].mat[13];
00268   cell[2] = Tpbc[0][1].mat[14];
00269   cell[3] = Tpbc[1][1].mat[12];
00270   cell[4] = Tpbc[1][1].mat[13];
00271   cell[5] = Tpbc[1][1].mat[14];
00272   cell[6] = Tpbc[2][1].mat[12];
00273   cell[7] = Tpbc[2][1].mat[13];
00274   cell[8] = Tpbc[2][1].mat[14];
00275 
00276   float len[3];
00277   len[0] = sqrtf(dot_prod(&cell[0], &cell[0]));
00278   len[1] = sqrtf(dot_prod(&cell[3], &cell[3]));
00279   len[2] = sqrtf(dot_prod(&cell[6], &cell[6]));
00280   //printf("len={%.3f %.3f %.3f}\n", len[0], len[1], len[2]);
00281 
00282   int i;
00283   float minlen = len[0];
00284   if (len[1] && len[1]<minlen) minlen = len[1];
00285   if (len[2] && len[2]<minlen) minlen = len[2];
00286   minlen--;
00287 
00288   // The algorithm works only for atoms in adjacent neighbor cells.
00289   if (!box && (cutoff[0]>=len[0] || cutoff[1]>=len[1] || cutoff[2]>=len[2])) {
00290     return MEASURE_ERR_BADCUTOFF;
00291   }
00292 
00293   bool bigrim = 1;
00294   float corecell[9];
00295   float diag[3];
00296   float origin[3];
00297   memset(origin, 0, 3L*sizeof(float));
00298   Matrix4 M_norm;
00299 
00300   if (box) {
00301     // Get the matrix M_norm that transforms all atoms inside the 
00302     // unit cell into the normalized unitcell spanned by 
00303     // {1/len[0] 0 0} {0 1/len[1] 0} {0 0 1/len[2]}.
00304     bigrim = 1;
00305 
00306     float vtmp[3];
00307     vec_add(vtmp, &cell[0], &cell[3]);
00308     vec_add(diag, &cell[6], vtmp);
00309     //printf("diag={%.3f %.3f %.3f}\n", diag[0], diag[1], diag[2]);
00310 
00311     // Finally we need to apply the translation of the cell origin
00312     vec_copy(origin, center);
00313     vec_scaled_add(origin, -0.5, &cell[0]);
00314     vec_scaled_add(origin, -0.5, &cell[3]);
00315     vec_scaled_add(origin, -0.5, &cell[6]);
00316     vec_negate(origin, origin);
00317     //printf("origin={%.3f %.3f %.3f}\n", origin[0], origin[1], origin[2]);
00318 
00319   } else if (2.0f*cutoff[0]<minlen && 2.0f*cutoff[1]<minlen && 2.0f*cutoff[2]<minlen) {
00320     // The cutoff must not be larger than half of the smallest cell dimension
00321     // otherwise we would have to use a less efficient algorithm.
00322 
00323     // Get the matrix M_norm that transforms all atoms inside the 
00324     // corecell into the orthonormal unitcell spanned by {1 0 0} {0 1 0} {0 0 1}.
00325     // The corecell ist the pbc cell minus cutoffs for each dimension.
00326     vec_scale(&corecell[0], (len[0]-cutoff[0])/len[0], &cell[0]);
00327     vec_scale(&corecell[3], (len[1]-cutoff[1])/len[1], &cell[3]);
00328     vec_scale(&corecell[6], (len[2]-cutoff[2])/len[2], &cell[6]);
00329     get_transform_to_orthonormal_cell(corecell, center, M_norm);
00330     //printf("Using algorithm for small PBC environment.\n");
00331 
00332   } else {
00333     // Get the matrix M_norm that transforms all atoms inside the 
00334     // unit cell into the orthonormal unitcell spanned by {1 0 0} {0 1 0} {0 0 1}.
00335     get_transform_to_orthonormal_cell(cell, center, M_norm);
00336 
00337     bigrim = 1;
00338     //printf("Using algorithm for large PBC environment.\n");
00339   }
00340 
00341   // In case the molecule was aligned our pbc cell is rotated and shifted.
00342   // In order to transform a point P into the orthonormal cell (P') it 
00343   // first has to be unaligned (the inverse of the alignment):
00344   // P' = M_norm * (alignment^-1) * P
00345   Matrix4 alignmentinv(*alignment);
00346   alignmentinv.inverse();
00347   Matrix4 M_coretransform(M_norm);
00348   M_coretransform.multmatrix(alignmentinv);
00349 
00350   //printf("alignment = \n");
00351   //print_Matrix4(alignment);
00352 
00353   // Similarly if we want to transform a point P into its image P' we
00354   // first have to unalign it, then apply the PBC translation and 
00355   // finally realign:
00356   // P' = alignment * Tpbc * (alignment^-1) * P
00357   //      `-------------v--------------'
00358   //                transform
00359   int j, u;
00360   Matrix4 Tpbc_aligned[3][2];
00361   if (!box) {
00362     for (i=0; i<3; i++) {
00363       for (j=0; j<2; j++) {
00364         Tpbc_aligned[i][j].loadmatrix(*alignment);
00365         Tpbc_aligned[i][j].multmatrix(Tpbc[i][j]);
00366         Tpbc_aligned[i][j].multmatrix(alignmentinv);
00367       }
00368     }
00369   }
00370 
00371   Matrix4 M[3];
00372   float *coords = ts->pos;
00373   float *coor;
00374   float orthcoor[3], wrapcoor[3];
00375 
00376   //printf("cutoff={%.3f %.3f %.3f}\n", cutoff[0], cutoff[1], cutoff[2]);
00377 
00378   if (box) {
00379     float min_coord[3], max_coord[3];
00380     // Increase box by cutoff
00381     vec_sub(min_coord, box,   cutoff);
00382     vec_add(max_coord, box+3, cutoff);
00383     //printf("Wrapping atoms into rectangular bounding box.\n");
00384     //printf("min_coord={%.3f %.3f %.3f}\n", min_coord[0], min_coord[1], min_coord[2]);
00385     //printf("max_coord={%.3f %.3f %.3f}\n", max_coord[0], max_coord[1], max_coord[2]);
00386     vec_add(min_coord, min_coord, origin);
00387     vec_add(max_coord, max_coord, origin);
00388 
00389     float testcoor[9];
00390     int idx, k;
00391     // Loop over all atoms
00392     for (idx=0; idx<ts->num; idx++) { 
00393       coor = coords+3L*idx;
00394 
00395       // Apply the inverse alignment transformation
00396       // to the current test point.
00397       M_coretransform.multpoint3d(coor, orthcoor);
00398 
00399       // Loop over all 26 neighbor cells
00400       // x
00401       for (i=-1; i<=1; i++) {
00402         // Choose the direction of translation
00403         if      (i>0) M[0].loadmatrix(Tpbc[0][1]);
00404         else if (i<0) M[0].loadmatrix(Tpbc[0][0]);
00405         else          M[0].identity();
00406         // Translate the unaligned atom
00407         M[0].multpoint3d(orthcoor, testcoor);
00408 
00409         // y
00410         for (j=-1; j<=1; j++) {
00411           // Choose the direction of translation
00412           if      (j>0) M[1].loadmatrix(Tpbc[1][1]);
00413           else if (j<0) M[1].loadmatrix(Tpbc[1][0]);
00414           else          M[1].identity();
00415           // Translate the unaligned atom
00416           M[1].multpoint3d(testcoor, testcoor+3);
00417 
00418           // z
00419           for (k=-1; k<=1; k++) {
00420             if(i==0 && j==0 && k==0) continue;
00421 
00422             // Choose the direction of translation
00423             if      (k>0) M[2].loadmatrix(Tpbc[2][1]);
00424             else if (k<0) M[2].loadmatrix(Tpbc[2][0]);
00425             else          M[2].identity();
00426             // Translate the unaligned atom
00427             M[2].multpoint3d(testcoor+3, testcoor+6);
00428 
00429             // Realign atom
00430             alignment->multpoint3d(testcoor+6, wrapcoor);
00431 
00432             vec_add(testcoor+6, wrapcoor, origin);
00433             if (testcoor[6]<min_coord[0] || testcoor[6]>max_coord[0]) continue;
00434             if (testcoor[7]<min_coord[1] || testcoor[7]>max_coord[1]) continue;
00435             if (testcoor[8]<min_coord[2] || testcoor[8]>max_coord[2]) continue;
00436 
00437             // Atom is inside cutoff, add it to the list            
00438             extcoord_array->append3(&wrapcoor[0]);
00439             indexmap_array->append(idx);
00440           }
00441         }
00442       }
00443     }
00444 
00445   } else if (bigrim) {
00446     // This is the more general but slower algorithm.
00447     // We loop over all atoms, move each atom to all 26 neighbor cells
00448     // and check if it lies inside cutoff
00449     float min_coord[3], max_coord[3];
00450     min_coord[0] = -cutoff[0]/len[0];
00451     min_coord[1] = -cutoff[1]/len[1];
00452     min_coord[2] = -cutoff[2]/len[2];
00453     max_coord[0] = 1.0f + cutoff[0]/len[0];
00454     max_coord[1] = 1.0f + cutoff[1]/len[1];
00455     max_coord[2] = 1.0f + cutoff[2]/len[2];
00456 
00457     float testcoor[3];
00458     int idx, k;
00459     // Loop over all atoms
00460     for (idx=0; idx<ts->num; idx++) { 
00461       coor = coords+3L*idx;
00462 
00463       // Apply the PBC --> orthonormal unitcell transformation
00464       // to the current test point.
00465       M_coretransform.multpoint3d(coor, orthcoor);
00466 
00467       // Loop over all 26 neighbor cells
00468       // x
00469       for (i=-1; i<=1; i++) {
00470         testcoor[0] = orthcoor[0]+(float)(i);
00471         if (testcoor[0]<min_coord[0] || testcoor[0]>max_coord[0]) continue;
00472 
00473         // Choose the direction of translation
00474         if      (i>0) M[0].loadmatrix(Tpbc_aligned[0][1]);
00475         else if (i<0) M[0].loadmatrix(Tpbc_aligned[0][0]);
00476         else          M[0].identity();
00477 
00478         // y
00479         for (j=-1; j<=1; j++) {
00480           testcoor[1] = orthcoor[1]+(float)(j);
00481           if (testcoor[1]<min_coord[1] || testcoor[1]>max_coord[1]) continue;
00482 
00483           // Choose the direction of translation
00484           if      (j>0) M[1].loadmatrix(Tpbc_aligned[1][1]);
00485           else if (j<0) M[1].loadmatrix(Tpbc_aligned[1][0]);
00486           else          M[1].identity();
00487 
00488           // z
00489           for (k=-1; k<=1; k++) {
00490             testcoor[2] = orthcoor[2]+(float)(k);
00491             if (testcoor[2]<min_coord[2] || testcoor[2]>max_coord[2]) continue;
00492 
00493             if(i==0 && j==0 && k==0) continue;
00494 
00495             // Choose the direction of translation
00496             if      (k>0) M[2].loadmatrix(Tpbc_aligned[2][1]);
00497             else if (k<0) M[2].loadmatrix(Tpbc_aligned[2][0]);
00498             else          M[2].identity();
00499 
00500             M[0].multpoint3d(coor, wrapcoor);
00501             M[1].multpoint3d(wrapcoor, wrapcoor);
00502             M[2].multpoint3d(wrapcoor, wrapcoor);
00503 
00504             // Atom is inside cutoff, add it to the list            
00505             extcoord_array->append3(&wrapcoor[0]);
00506             indexmap_array->append(idx);
00507           }
00508         }
00509       }
00510     }
00511   
00512   } else {
00513     Matrix4 Mtmp;
00514 
00515     for (i=0; i < ts->num; i++) { 
00516       // Apply the PBC --> orthonormal unitcell transformation
00517       // to the current test point.
00518       M_coretransform.multpoint3d(coords+3L*i, orthcoor);
00519 
00520       // Determine in which cell we are.
00521       int cellindex[3];    
00522       if      (orthcoor[0]<0) cellindex[0] = -1;
00523       else if (orthcoor[0]>1) cellindex[0] =  1;
00524       else                    cellindex[0] =  0;
00525       if      (orthcoor[1]<0) cellindex[1] = -1;
00526       else if (orthcoor[1]>1) cellindex[1] =  1;
00527       else                    cellindex[1] =  0;
00528       if      (orthcoor[2]<0) cellindex[2] = -1;
00529       else if (orthcoor[2]>1) cellindex[2] =  1;
00530       else                    cellindex[2] =  0;
00531 
00532       // All zero means we're inside the core --> no image.
00533       if (!cellindex[0] && !cellindex[1] && !cellindex[2]) continue;
00534 
00535       // Choose the direction of translation
00536       if      (orthcoor[0]<0) M[0].loadmatrix(Tpbc_aligned[0][1]);
00537       else if (orthcoor[0]>1) M[0].loadmatrix(Tpbc_aligned[0][0]);
00538       if      (orthcoor[1]<0) M[1].loadmatrix(Tpbc_aligned[1][1]);
00539       else if (orthcoor[1]>1) M[1].loadmatrix(Tpbc_aligned[1][0]);
00540       if      (orthcoor[2]<0) M[2].loadmatrix(Tpbc_aligned[2][1]);
00541       else if (orthcoor[2]>1) M[2].loadmatrix(Tpbc_aligned[2][0]);
00542 
00543       // Create wrapped copies of the atom:
00544       // x, y, z planes
00545       coor = coords+3L*i;
00546       for (u=0; u<3; u++) {
00547         if (cellindex[u] && cutoff[u]) {
00548           M[u].multpoint3d(coor, wrapcoor);
00549           extcoord_array->append3(&wrapcoor[0]);
00550           indexmap_array->append(i);
00551         }
00552       }
00553 
00554       Mtmp = M[0];
00555 
00556       // xy edge
00557       if (cellindex[0] && cellindex[1] && cutoff[0] && cutoff[1]) {
00558         M[0].multmatrix(M[1]);
00559         M[0].multpoint3d(coor, wrapcoor);
00560         extcoord_array->append3(&wrapcoor[0]);
00561         indexmap_array->append(i);
00562       }
00563 
00564       // yz edge
00565       if (cellindex[1] && cellindex[2] && cutoff[1] && cutoff[2]) {
00566         M[1].multmatrix(M[2]);
00567         M[1].multpoint3d(coor, wrapcoor);
00568         extcoord_array->append3(&wrapcoor[0]);
00569         indexmap_array->append(i);
00570       }
00571 
00572       // zx edge
00573       if (cellindex[0] && cellindex[2] && cutoff[0] && cutoff[2]) {
00574         M[2].multmatrix(Mtmp);
00575         M[2].multpoint3d(coor, wrapcoor);
00576         extcoord_array->append3(&wrapcoor[0]);
00577         indexmap_array->append(i);
00578       }
00579 
00580       // xyz corner
00581       if (cellindex[0] && cellindex[1] && cellindex[2]) {
00582         M[1].multmatrix(Mtmp);
00583         M[1].multpoint3d(coor, wrapcoor);
00584         extcoord_array->append3(&wrapcoor[0]);
00585         indexmap_array->append(i);
00586       }
00587 
00588     }
00589 
00590   } // endif
00591 
00592   // If a selection was provided we select extcoords
00593   // within cutoff of the original selection:
00594   if (sel) {
00595     int numext = sel->selected+indexmap_array->num();
00596     float *extcoords = new float[3L*numext];
00597     int   *indexmap  = new int[numext];
00598     int   *others    = new int[numext];
00599     memset(others, 0, numext);
00600 
00601     // Use the largest given cutoff
00602     float maxcutoff = cutoff[0];
00603     for (i=1; i<3; i++) {
00604       if (cutoff[i]>maxcutoff) maxcutoff = cutoff[i];
00605     }
00606 
00607     // Prepare C-array of coordinates for find_within()
00608     j=0;
00609     for (i=0; i < sel->num_atoms; i++) { 
00610       if (!sel->on[i]) continue; //atom is not selected
00611       extcoords[3L*j]   = coords[3L*i];
00612       extcoords[3L*j+1] = coords[3L*i+1];
00613       extcoords[3L*j+2] = coords[3L*i+2];
00614       indexmap[j] = i;
00615       others[j++] = 1;
00616     }
00617     for (i=0; i<indexmap_array->num(); i++) {
00618       extcoords[3L*j]   = (*extcoord_array)[3L*i];
00619       extcoords[3L*j+1] = (*extcoord_array)[3L*i+1];
00620       extcoords[3L*j+2] = (*extcoord_array)[3L*i+2];
00621       indexmap[j] = (*indexmap_array)[i];
00622       others[j++] = 0;
00623     }
00624 
00625     // Initialize flags array to true, find_within() results are AND'd/OR'd in.
00626     int *flgs   = new int[numext];
00627     for (i=0; i<numext; i++) {
00628       flgs[i] = 1;
00629     }
00630 
00631     // Find coordinates from extcoords that are within cutoff of the ones
00632     // with flagged in 'others' and set the flgs accordingly:
00633     find_within(extcoords, flgs, others, numext, maxcutoff);
00634 
00635     extcoord_array->clear();
00636     indexmap_array->clear();
00637     for (i=sel->selected; i<numext; i++) {
00638       if (!flgs[i]) continue;
00639 
00640       extcoord_array->append3(&extcoords[3L*i]);
00641       indexmap_array->append(indexmap[i]);
00642     }
00643 
00644   }
00645 
00646   return MEASURE_NOERR;
00647 }  
00648 
00649 // Computes the rectangular bounding box for the PBC cell.
00650 // If the molecule was rotated/moved you can supply the transformation
00651 // matrix and you'll get the bounding box of the transformed cell.
00652 int compute_pbcminmax(MoleculeList *mlist, int molid, int frame, 
00653                const float *center, const Matrix4 *transform,
00654                float *min, float *max) {
00655   Molecule *mol = mlist->mol_from_id(molid);
00656   if( !mol )
00657     return MEASURE_ERR_NOMOLECULE;
00658 
00659   Timestep *ts = mol->get_frame(frame);
00660   if (!ts) return MEASURE_ERR_NOFRAMES;
00661 
00662   // Get the displacement vectors (in form of translation matrices)
00663   Matrix4 Tpbc[3];
00664   ts->get_transforms(Tpbc[0], Tpbc[1], Tpbc[2]);
00665 
00666   // Construct the cell spanning vectors
00667   float cell[9];
00668   cell[0] = Tpbc[0].mat[12];
00669   cell[1] = Tpbc[0].mat[13];
00670   cell[2] = Tpbc[0].mat[14];
00671   cell[3] = Tpbc[1].mat[12];
00672   cell[4] = Tpbc[1].mat[13];
00673   cell[5] = Tpbc[1].mat[14];
00674   cell[6] = Tpbc[2].mat[12];
00675   cell[7] = Tpbc[2].mat[13];
00676   cell[8] = Tpbc[2].mat[14];
00677 
00678 #if 0
00679   float len[3];
00680   len[0] = sqrtf(dot_prod(&cell[0], &cell[0]));
00681   len[1] = sqrtf(dot_prod(&cell[3], &cell[3]));
00682   len[2] = sqrtf(dot_prod(&cell[6], &cell[6]));
00683 #endif
00684 
00685   // Construct all 8 corners (nodes) of the bounding box
00686   float node[8*3];
00687   int n=0;
00688   float i, j, k;
00689   for (i=-0.5; i<1.f; i+=1.f) {
00690     for (j=-0.5; j<1.f; j+=1.f) {
00691       for (k=-0.5; k<1.f; k+=1.f) {
00692         // Apply the translation of the origin
00693         vec_copy(node+3L*n, center);
00694         vec_scaled_add(node+3L*n, i, &cell[0]);
00695         vec_scaled_add(node+3L*n, j, &cell[3]);
00696         vec_scaled_add(node+3L*n, k, &cell[6]);
00697 
00698         // Apply global alignment transformation
00699         transform->multpoint3d(node+3L*n, node+3L*n);
00700         n++;
00701       }
00702     }
00703   }
00704 
00705   // Find minmax coordinates of all corners
00706   for (n=0; n<8; n++) {
00707     if (!n || node[3L*n  ]<min[0])  min[0] = node[3L*n];
00708     if (!n || node[3L*n+1]<min[1])  min[1] = node[3L*n+1];
00709     if (!n || node[3L*n+2]<min[2])  min[2] = node[3L*n+2];
00710     if (!n || node[3L*n  ]>max[0])  max[0] = node[3L*n];
00711     if (!n || node[3L*n+1]>max[1])  max[1] = node[3L*n+1];
00712     if (!n || node[3L*n+2]>max[2])  max[2] = node[3L*n+2];
00713   }
00714 
00715   return MEASURE_NOERR;
00716 }

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