00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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,
00051
00052 float* xyz,
00053
00054 int natoms2,
00055
00056 float* xyz2,
00057
00058 float* cell,
00059 float* hist,
00060
00061 int maxbin,
00062 float rmin,
00063 float delr)
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;
00142 int duplicates=0;
00143 float rmin = 0.0f;
00144
00145
00146 a=b=c=9999999.0;
00147 alpha=beta=gamma=90.0;
00148
00149
00150 framecntr[0]=framecntr[1]=framecntr[2]=0;
00151
00152
00153
00154 if (!sel1 || !sel2 ) {
00155 return MEASURE_ERR_NOSEL;
00156 }
00157
00158
00159
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
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
00185 ts = sel1->timestep(mlist);
00186 frame=last;
00187 } else {
00188 ts = mymol->get_frame(frame);
00189 }
00190
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
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
00205 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0)) {
00206 isortho=0;
00207 }
00208 }
00209 } else {
00210
00211 isortho=1;
00212 a=b=c=9999999.0;
00213 alpha=beta=gamma=90.0;
00214 }
00215
00216
00217 if (!isortho) {
00218 msgErr << "measure rdf: only orthorhombic cells are supported (for now)." << sendmsg;
00219 return MEASURE_ERR_GENERAL;
00220 }
00221
00222
00223 for (i=0; i<count_h; ++i) {
00224 gofr[i] = numint[i] = histog[i] = 0.0;
00225 }
00226
00227
00228
00229
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
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);
00246 }
00247
00248 if (usepbc) {
00249
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
00259 float boxby2[3];
00260 boxby2[0] = 0.5f * a;
00261 boxby2[1] = 0.5f * b;
00262 boxby2[2] = 0.5f * c;
00263
00264
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
00273
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
00283
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
00306 pbccell[0]=a;
00307 pbccell[1]=b;
00308 pbccell[2]=c;
00309
00310
00311 memset(lhist, 0, count_h * sizeof(float));
00312
00313 if (isortho && sel1->selected && sel2->selected) {
00314
00315
00316 int rc=-1;
00317 #if defined(VMDCUDA)
00318 if (!getenv("VMDNOCUDA") && (app->cuda != NULL)) {
00319
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
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];
00343 } else {
00344 ++framecntr[1];
00345 }
00346 ++framecntr[0];
00347
00348 #if 0
00349
00350
00351
00352
00353
00354
00355 lhist[0] -= duplicates;
00356 #endif
00357
00358
00359
00360
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
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 {
00384 pair_dens = 30.0 * (double)sel1->selected /
00385 ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
00386 }
00387 }
00388
00389
00390 for (i=0; i<h_max; ++i) {
00391
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
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