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