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 }