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

MeasureRDF.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: MeasureRDF.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.17 $       $Date: 2019/01/17 21:21:00 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code to compute radial distribution functions for MD trajectories.
00019  *
00020  ***************************************************************************/
00021 
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include "Measure.h"
00026 #include "AtomSel.h"
00027 #include "utilities.h"
00028 #include "ResizeArray.h"
00029 #include "MoleculeList.h"
00030 #include "Inform.h"
00031 #include "Timestep.h"
00032 #include "VMDApp.h"
00033 #include "WKFThreads.h"
00034 #include "WKFUtils.h"
00035 #include "CUDAAccel.h"
00036 #include "CUDAKernels.h"
00037 
00038 #define MIN(X,Y) (((X)<(Y))? (X) : (Y))
00039 #define MAX(X,Y) (((X)>(Y))? (X) : (Y))
00040 
00045 static inline double spherical_cap(const double &radius, const double &boxby2) {
00046   return (VMD_PI / 3.0 * (radius - boxby2) * (radius - boxby2)
00047           * ( 2.0 * radius + boxby2));
00048 }
00049 
00050 void rdf_cpu(int natoms1,     // array of the number of atoms in
00051                               // selection 1 in each frame. 
00052              float* xyz,      // coordinates of first selection.
00053                               // [natoms1][3]
00054              int natoms2,     // array of the number of atoms in
00055                               // selection 2 in each frame. 
00056              float* xyz2,     // coordinates of selection 2.
00057                               // [natoms2][3]
00058              float* cell,     // the cell x y and z dimensions [3]
00059              float* hist,     // the histograms, 1 per block
00060                               // [ncudablocks][maxbin]
00061              int maxbin,      // the number of bins in the histogram
00062              float rmin,      // the minimum value of the first bin
00063              float delr)      // the width of each bin
00064 {
00065   int iatom, jatom, ibin;
00066   float rij, rxij, rxij2, x1, y1, z1, x2, y2, z2;
00067   float cellx, celly, cellz;
00068   int *ihist = new int[maxbin];
00069 
00070   for (ibin=0; ibin<maxbin; ibin++) {
00071     ihist[ibin]=0;
00072   }
00073 
00074   cellx = cell[0];
00075   celly = cell[1];
00076   cellz = cell[2];
00077 
00078   for (iatom=0; iatom<natoms1; iatom++) {
00079     long addr = 3L * iatom;
00080     xyz[addr    ] = fmodf(xyz[addr    ], cellx);
00081     xyz[addr + 1] = fmodf(xyz[addr + 1], celly);
00082     xyz[addr + 2] = fmodf(xyz[addr + 2], cellz);
00083   }
00084 
00085   for (iatom=0; iatom<natoms2; iatom++) {
00086     long addr = 3L * iatom;
00087     xyz2[addr    ] = fmodf(xyz2[addr    ], cellx);
00088     xyz2[addr + 1] = fmodf(xyz2[addr + 1], celly);
00089     xyz2[addr + 2] = fmodf(xyz2[addr + 2], cellz);
00090   }
00091 
00092   for (iatom=0; iatom<natoms1; iatom++) {
00093     x1 = xyz[3L * iatom    ];
00094     y1 = xyz[3L * iatom + 1];
00095     z1 = xyz[3L * iatom + 2];
00096     for (jatom=0;jatom<natoms2;jatom++) {
00097       x2 = xyz2[3L * jatom    ];
00098       y2 = xyz2[3L * jatom + 1];
00099       z2 = xyz2[3L * jatom + 2];
00100 
00101       rxij = fabsf(x1 - x2);
00102       rxij2 = cellx - rxij;
00103       rxij = MIN(rxij, rxij2);
00104       rij = rxij * rxij;
00105 
00106       rxij = fabsf(y1 - y2);
00107       rxij2 = celly - rxij;
00108       rxij = MIN(rxij, rxij2);
00109       rij += rxij * rxij;
00110 
00111       rxij = fabsf(z1 - z2);
00112       rxij2 = cellz - rxij;
00113       rxij = MIN(rxij, rxij2);
00114       rij += rxij * rxij;
00115 
00116       rij = sqrtf(rij);
00117 
00118       ibin = (int)floorf((rij-rmin)/delr);
00119       if (ibin<maxbin && ibin>=0) {
00120         ++ihist[ibin];
00121       }
00122     }
00123   }
00124 
00125   delete [] ihist;
00126 }
00127 
00128 
00129 
00130 
00131 int measure_rdf(VMDApp *app,
00132                 AtomSel *sel1, AtomSel *sel2, MoleculeList *mlist,
00133                 const int count_h, double *gofr, 
00134                 double *numint, double *histog,
00135                 const float delta, int first, int last, int step, 
00136                 int *framecntr, int usepbc, int selupdate) {
00137   int i, j, frame;
00138 
00139   float a, b, c, alpha, beta, gamma;
00140   float pbccell[3];
00141   int isortho=0;     // orthogonal unit cell not assumed by default.
00142   int duplicates=0;  // counter for duplicate atoms in both selections.
00143   float rmin = 0.0f; // min distance to histogram
00144 
00145   // initialize a/b/c/alpha/beta/gamma to arbitrary defaults to please the compiler.
00146   a=b=c=9999999.0;
00147   alpha=beta=gamma=90.0;
00148 
00149   // reset counter for total, skipped, and _orth processed frames.
00150   framecntr[0]=framecntr[1]=framecntr[2]=0;
00151 
00152   // First round of sanity checks.
00153   // neither list can be undefined
00154   if (!sel1 || !sel2 ) {
00155     return MEASURE_ERR_NOSEL;
00156   }
00157 
00158   // make sure that both selections are from the same molecule
00159   // so that we know that PBC unit cell info is the same for both
00160   if (sel2->molid() != sel1->molid()) {
00161     return MEASURE_ERR_MISMATCHEDMOLS;
00162   }
00163 
00164   Molecule *mymol = mlist->mol_from_id(sel1->molid());
00165   int maxframe = mymol->numframes() - 1;
00166   int nframes = 0;
00167 
00168   if (last == -1)
00169     last = maxframe;
00170 
00171   if ((last < first) || (last < 0) || (step <=0) || (first < -1)
00172       || (maxframe < 0) || (last > maxframe)) {
00173       msgErr << "measure rdf: bad frame range given."
00174              << " max. allowed frame#: " << maxframe << sendmsg;
00175     return MEASURE_ERR_BADFRAMERANGE;
00176   }
00177 
00178   // test for non-orthogonal PBC cells, zero volume cells, etc.
00179   if (usepbc) {
00180     for (isortho=1, nframes=0, frame=first; frame <=last; ++nframes, frame += step) {
00181       const Timestep *ts;
00182 
00183       if (first == -1) {
00184         // use current frame only. don't loop.
00185         ts = sel1->timestep(mlist);
00186         frame=last;
00187       } else {
00188         ts = mymol->get_frame(frame);
00189       }
00190       // get periodic cell information for current frame
00191       a = ts->a_length;
00192       b = ts->b_length;
00193       c = ts->c_length;
00194       alpha = ts->alpha;
00195       beta = ts->beta;
00196       gamma = ts->gamma;
00197 
00198       // check validity of PBC cell side lengths
00199       if (fabsf(a*b*c) < 0.0001) {
00200         msgErr << "measure rdf: unit cell volume is zero." << sendmsg;
00201         return MEASURE_ERR_GENERAL;
00202       }
00203 
00204       // check PBC unit cell shape to select proper low level algorithm.
00205       if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0)) {
00206         isortho=0;
00207       }
00208     }
00209   } else {
00210     // initialize a/b/c/alpha/beta/gamma to arbitrary defaults
00211     isortho=1;
00212     a=b=c=9999999.0;
00213     alpha=beta=gamma=90.0;
00214   }
00215 
00216   // until we can handle non-orthogonal periodic cells, this is fatal
00217   if (!isortho) {
00218     msgErr << "measure rdf: only orthorhombic cells are supported (for now)." << sendmsg;
00219     return MEASURE_ERR_GENERAL;
00220   }
00221 
00222   // clear the result arrays
00223   for (i=0; i<count_h; ++i) {
00224     gofr[i] = numint[i] = histog[i] = 0.0;
00225   }
00226 
00227   // pre-allocate coordinate buffers of the max size we'll
00228   // ever need, so we don't have to reallocate if/when atom
00229   // selections are updated on-the-fly
00230   float *sel1coords = new float[3L*sel1->num_atoms];
00231   float *sel2coords = new float[3L*sel2->num_atoms];
00232   float *lhist = new float[count_h];
00233 
00234   for (nframes=0,frame=first; frame <=last; ++nframes, frame += step) {
00235     const Timestep *ts1, *ts2;
00236 
00237     if (frame  == -1) {
00238       // use current frame only. don't loop.
00239       ts1 = sel1->timestep(mlist);
00240       ts2 = sel2->timestep(mlist);
00241       frame=last;
00242     } else {
00243       sel1->which_frame = frame;
00244       sel2->which_frame = frame;
00245       ts1 = ts2 = mymol->get_frame(frame); // requires sels from same mol
00246     }
00247 
00248     if (usepbc) {
00249       // get periodic cell information for current frame
00250       a     = ts1->a_length;
00251       b     = ts1->b_length;
00252       c     = ts1->c_length;
00253       alpha = ts1->alpha;
00254       beta  = ts1->beta;
00255       gamma = ts1->gamma;
00256     }
00257 
00258     // compute half periodic cell size
00259     float boxby2[3];
00260     boxby2[0] = 0.5f * a;
00261     boxby2[1] = 0.5f * b;
00262     boxby2[2] = 0.5f * c;
00263 
00264     // update the selections if the user desires it
00265     if (selupdate) {
00266       if (sel1->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
00267         msgErr << "measure rdf: failed to evaluate atom selection update";
00268       if (sel2->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
00269         msgErr << "measure rdf: failed to evaluate atom selection update";
00270     }
00271 
00272     // check for duplicate atoms in the two lists, as these will have
00273     // to be subtracted back out of the first histogram slot
00274     if (sel2->molid() == sel1->molid()) {
00275       int i;
00276       for (i=0, duplicates=0; i<sel1->num_atoms; ++i) {
00277         if (sel1->on[i] && sel2->on[i])
00278           ++duplicates;
00279       }
00280     }
00281 
00282     // copy selected atoms to the two coordinate lists
00283     // requires that selections come from the same molecule
00284     const float *framepos = ts1->pos;
00285     for (i=0, j=0; i<sel1->num_atoms; ++i) {
00286       if (sel1->on[i]) {
00287         long a = i*3L;
00288         sel1coords[j    ] = framepos[a    ];
00289         sel1coords[j + 1] = framepos[a + 1];
00290         sel1coords[j + 2] = framepos[a + 2];
00291         j+=3;
00292       }
00293     }
00294     framepos = ts2->pos;
00295     for (i=0, j=0; i<sel2->num_atoms; ++i) {
00296       if (sel2->on[i]) {
00297         long a = i*3L;
00298         sel2coords[j    ] = framepos[a    ];
00299         sel2coords[j + 1] = framepos[a + 1];
00300         sel2coords[j + 2] = framepos[a + 2];
00301         j+=3;
00302       }
00303     }
00304 
00305     // copy unit cell information
00306     pbccell[0]=a;
00307     pbccell[1]=b;
00308     pbccell[2]=c;
00309 
00310     // clear the histogram for this frame
00311     memset(lhist, 0, count_h * sizeof(float));
00312 
00313     if (isortho && sel1->selected && sel2->selected) {
00314       // do the rdf calculation for orthogonal boxes.
00315       // XXX. non-orthogonal box not supported yet. detected and handled above.
00316       int rc=-1;
00317 #if defined(VMDCUDA)
00318       if (!getenv("VMDNOCUDA") && (app->cuda != NULL)) {
00319 //        msgInfo << "Running multi-GPU RDF calc..." << sendmsg;
00320         rc=rdf_gpu(app->cuda->get_cuda_devpool(),
00321                    usepbc,
00322                    sel1->selected, sel1coords,
00323                    sel2->selected, sel2coords, 
00324                    pbccell,
00325                    lhist,
00326                    count_h,
00327                    rmin,
00328                    delta);
00329       } 
00330 #endif
00331       if (rc != 0) {
00332 //        msgInfo << "Running single-core CPU RDF calc..." << sendmsg;
00333         rdf_cpu(sel1->selected, sel1coords,
00334                 sel2->selected, sel2coords, 
00335                 pbccell,
00336                 lhist,
00337                 count_h,
00338                 rmin,
00339                 delta);
00340       }
00341 
00342       ++framecntr[2]; // frame processed with rdf algorithm
00343     } else {
00344       ++framecntr[1]; // frame skipped
00345     }
00346     ++framecntr[0];   // total frames.
00347 
00348 #if 0
00349     // XXX elimination of duplicates is now handled within the 
00350     //     GPU kernels themselves, so we do not need to subtract them
00351     //     off during the histogram normalization calculations.
00352     // Correct the first histogram slot for the number of atoms that are
00353     // present in both lists. they'll end up in the first histogram bin.
00354     // we subtract only from the first thread histogram which is always defined.
00355     lhist[0] -= duplicates;
00356 #endif
00357 
00358     // in case of going 'into the edges', we should cut
00359     // off the part that is not properly normalized to
00360     // not confuse people that don't know about this.
00361     int h_max=count_h;
00362     float smallside=a;
00363     if (isortho && usepbc) {
00364       if(b < smallside) {
00365         smallside=b;
00366       }
00367       if(c < smallside) {
00368         smallside=c;
00369       }
00370       h_max=(int) (sqrtf(0.5f)*smallside/delta) +1;
00371       if (h_max > count_h) {
00372         h_max=count_h;
00373       }
00374     }
00375 
00376     // compute normalization function.
00377     double all=0.0;
00378     double pair_dens = 0.0;
00379     
00380     if (sel1->selected && sel2->selected) {
00381       if (usepbc) {
00382         pair_dens = a * b * c / ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
00383       } else { // assume a particle volume of 30 \AA^3 (~ 1 water).
00384         pair_dens = 30.0 * (double)sel1->selected /
00385           ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
00386       }
00387     }
00388 
00389     // XXX for orthogonal boxes, we can reduce this to rmax < sqrt(0.5)*smallest side
00390     for (i=0; i<h_max; ++i) {
00391       // radius of inner and outer sphere that form the spherical slice
00392       double r_in  = delta * (double)i;
00393       double r_out = delta * (double)(i+1);
00394       double slice_vol = 4.0 / 3.0 * VMD_PI
00395         * ((r_out * r_out * r_out) - (r_in * r_in * r_in));
00396 
00397       if (isortho && usepbc) {
00398         // add correction for 0.5*box < r <= sqrt(0.5)*box
00399         if (r_out > boxby2[0]) {
00400           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[0]);
00401         }
00402         if (r_out > boxby2[1]) {
00403           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[1]);
00404         }
00405         if (r_out > boxby2[2]) {
00406           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[2]);
00407         }
00408         if (r_in > boxby2[0]) {
00409           slice_vol += 2.0 * spherical_cap(r_in, boxby2[0]);
00410         }
00411         if (r_in > boxby2[1]) {
00412           slice_vol += 2.0 * spherical_cap(r_in, boxby2[1]);
00413         }
00414         if (r_in > boxby2[2]) {
00415           slice_vol += 2.0 * spherical_cap(r_in, boxby2[2]);
00416         }
00417       }
00418 
00419       double normf = pair_dens / slice_vol;
00420       double histv = (double) lhist[i];
00421       gofr[i] += normf * histv;
00422       all     += histv;
00423       if (sel1->selected) {
00424         numint[i] += all / (double)(sel1->selected);
00425       }
00426       histog[i] += histv;
00427     }
00428   }
00429   delete [] sel1coords;
00430   delete [] sel2coords;
00431   delete [] lhist;
00432 
00433   double norm = 1.0 / (double) nframes;
00434   for (i=0; i<count_h; ++i) {
00435     gofr[i]   *= norm;
00436     numint[i] *= norm;
00437     histog[i] *= norm;
00438   }
00439 
00440   return MEASURE_NOERR;
00441 }
00442 

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