00001
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #define BENCH_PERIOD 1000
00023
00024
00025 #define GBIS_DEDR_FORCE 1
00026 #define GBIS_DEDA_FORCE 1
00027 #define GBIS_COUL_FORCE 1
00028
00029 #include "Vector.h"
00030 #include <limits>
00031 #include "InfoStream.h"
00032 #include "Node.h"
00033 #include "PatchMap.h"
00034 #include "PatchMap.inl"
00035 #include "AtomMap.h"
00036 #include "ComputeGBISser.h"
00037 #include "ComputeGBISserMgr.decl.h"
00038 #include "PatchMgr.h"
00039 #include "Molecule.h"
00040 #include "ReductionMgr.h"
00041 #include "ComputeMgr.h"
00042 #include "ComputeMgr.decl.h"
00043 #include "ComputeGBIS.inl"
00044
00045 #define MIN_DEBUG_LEVEL 3
00046 #include "Debug.h"
00047 #include "SimParameters.h"
00048 #include "WorkDistrib.h"
00049 #include "varsizemsg.h"
00050 #include <stdlib.h>
00051 #include <stdio.h>
00052 #include <errno.h>
00053 #include <time.h>
00054
00055 struct vect {
00056 BigReal x, y, z;
00057 };
00058 #ifndef COUL_CONST
00059 #define COUL_CONST 332.0636 // ke [kcal*Ang/e^2]
00060 #endif
00061
00062 #ifndef TA
00063 #define TA 0.333333333333333 // 1/3
00064 #define TB 0.4 // 2/5
00065 #define TC 0.428571428571428 // 3/7
00066 #define TD 0.444444444444444 // 4/9
00067 #define TE 0.454545454545454 // 5/11
00068 #define DA 1.333333333333333 // 4* 1/3
00069 #define DB 2.4 // 6* 2/5
00070 #define DC 3.428571428571428 // 8* 3/7
00071 #define DD 4.444444444444444 // 10*4/9
00072 #define DE 5.454545454545454 // 12*5/11
00073 #endif
00074
00075 inline void Phase2_PairSer(
00076
00077
00078 BigReal r,
00079 BigReal r2,
00080 BigReal r_i,
00081 BigReal qiqj,
00082 BigReal ai,
00083 BigReal aj,
00084 BigReal epsilon_p_i,
00085 BigReal epsilon_s_i,
00086 BigReal kappa,
00087 int exclij,
00088 BigReal scale14,
00089 int stat,
00090
00091
00092 BigReal & coulEij,
00093 BigReal & ddrCoulEij,
00094 BigReal & gbEij,
00095 BigReal & ddrGbEij,
00096 BigReal & dEdai,
00097 BigReal & dEdaj
00098 );
00099
00100 inline void CalcHPairSer (
00101 BigReal r,
00102 BigReal r2,
00103 BigReal ri,
00104 BigReal rc,
00105 BigReal ri0,
00106 BigReal rjs,
00107 BigReal rj0,
00108 BigReal ris,
00109 int & dij,
00110 int & dji,
00111 BigReal & dhij,
00112 BigReal & dhji);
00113 inline void CalcDHPairSer (
00114 BigReal r,
00115 BigReal r2,
00116 BigReal ri,
00117 BigReal rc,
00118 BigReal ri0,
00119 BigReal rjs,
00120 BigReal rj0,
00121 BigReal ris,
00122 int & dij,
00123 int & dji,
00124 BigReal & dhij,
00125 BigReal & dhji);
00126
00127
00128 struct ComputeGBISAtom {
00129 Position position;
00130 Real charge;
00131 Real rho;
00132 Real rho0;
00133 Real rhos;
00134 Mass mass;
00135
00136 int id;
00137 int vdwType;
00138 };
00139
00140 class GBISCoordMsg : public CMessage_GBISCoordMsg {
00141 public:
00142 int sourceNode;
00143 int numAtoms;
00144 int doSlow;
00145 ComputeGBISAtom *coord;
00146 int sequence;
00147 };
00148
00149 class GBISForceMsg : public CMessage_GBISForceMsg {
00150 public:
00151 BigReal gbInterEnergy;
00152 BigReal gbSelfEnergy;
00153 BigReal coulEnergy;
00154 ExtForce *force;
00155 ExtForce *slowForce;
00156 };
00157
00158 class ComputeGBISserMgr : public CBase_ComputeGBISserMgr {
00159 public:
00160 ComputeGBISserMgr();
00161 ~ComputeGBISserMgr();
00162
00163 void setCompute(ComputeGBISser *c) { gbisCompute = c; }
00164 void recvCoord(GBISCoordMsg *);
00165 void recvForce(GBISForceMsg *);
00166
00167 private:
00168 CProxy_ComputeGBISserMgr gbisProxy;
00169 ComputeGBISser *gbisCompute;
00170
00171 int numSources;
00172 int numArrived;
00173 GBISCoordMsg **coordMsgs;
00174 int numAtoms;
00175 ComputeGBISAtom *coord;
00176 ExtForce *force;
00177 ExtForce *slowForce;
00178 GBISForceMsg *oldmsg;
00179
00180 BigReal gbSelfEnergy;
00181 BigReal gbInterEnergy;
00182 BigReal coulEnergy;
00183 int timestep;
00184 clock_t t_start;
00185 clock_t t_stop;
00186 clock_t t1, t2;
00187 double totalnamdtime;
00188 double totalgbistime;
00189 double inittime;
00190 double psitime;
00191 double alphatime;
00192 double dEdrtime;
00193 double dEdasumtime;
00194 double dEdaprefixtime;
00195 double dEdalooptime;
00196 void calcGBISReg(int stat);
00197 int blockSize;
00198 int numBlocks;
00199 int loop1Iter, loop2Iter, loop3Iter;
00200 bool all2all;
00201 BigReal MinErr;
00202 BigReal MaxErr;
00203 int numBins;
00204 BigReal base;
00205 };
00206
00207 ComputeGBISserMgr::ComputeGBISserMgr() :
00208 gbisProxy(thisgroup), gbisCompute(0), numSources(0), numArrived(0),
00209 inittime(0),
00210 psitime(0),
00211 alphatime(0),
00212 dEdrtime(0),
00213 dEdasumtime(0),
00214 dEdaprefixtime(0),
00215 dEdalooptime(0), slowForce(0),
00216 coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0), timestep(0) {
00217 CkpvAccess(BOCclass_group).computeGBISserMgr = thisgroup;
00218 t_start = clock();
00219 t_stop = clock();
00220 all2all = false;
00221 }
00222
00223 ComputeGBISserMgr::~ComputeGBISserMgr() {
00224 for ( int i=0; i<numSources; ++i ) { delete coordMsgs[i]; }
00225 delete [] coordMsgs;
00226 delete [] coord;
00227 delete [] force;
00228 delete [] slowForce;
00229 delete oldmsg;
00230 }
00231
00232 ComputeGBISser::ComputeGBISser(ComputeID c) :
00233 ComputeHomePatches(c)
00234 {
00235 CProxy_ComputeGBISserMgr::ckLocalBranch(
00236 CkpvAccess(BOCclass_group).computeGBISserMgr)->setCompute(this);
00237
00238 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00239
00240 }
00241
00242 ComputeGBISser::~ComputeGBISser()
00243 {
00244 }
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 inline void CalcScaleSer (
00258 BigReal r,
00259 BigReal t,
00260 BigReal c,
00261 BigReal & s,
00262 BigReal & d ) {
00263
00264 if (r <= t) {
00265 s = 1;
00266 d = 0;
00267 } else if (r < c) {
00268
00269 BigReal ct = (c-t);
00270 BigReal ct_i = 1/ct;
00271 BigReal ct_i2 = ct_i*ct_i;
00272
00273
00274 BigReal rt = r - t;
00275 BigReal rt2 = rt*rt;
00276 BigReal omrtct2 = 1-rt2*ct_i2;
00277 s = omrtct2*omrtct2;
00278
00279 BigReal rc = r - c;
00280 d=s*4*rt/(rc*(rt+ct));
00281
00282 } else {
00283 s = 0;
00284 d = 0;
00285 }
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 void ComputeGBISserMgr::calcGBISReg(int stat) {
00297 t1 = clock();
00298
00299
00300
00301
00302
00303 Molecule *molecule = Node::Object()->molecule;
00304 SimParameters *simParams = Node::Object()->simParameters;
00305 BigReal epsilon_s = simParams->solvent_dielectric;
00306 BigReal epsilon_p = simParams->dielectric;
00307 BigReal epsilon_s_i = 1/epsilon_s;
00308 BigReal epsilon_p_i = 1/epsilon_p;
00309 BigReal rho_0 = simParams->coulomb_radius_offset;
00310 BigReal cutoff = simParams->cutoff;
00311 BigReal trans = simParams->switchingDist;
00312 BigReal fsMax = 1.728;
00313 BigReal a_cut = simParams->alpha_cutoff-fsMax;
00314
00315 BigReal a_cut2 = a_cut*a_cut;
00316 BigReal a_cut_ps = a_cut + fsMax;
00317 BigReal cutoff2 = cutoff*cutoff;
00318 BigReal a_cut_ps2 = a_cut_ps*a_cut_ps;
00319 BigReal delta = simParams->gbis_delta;
00320 BigReal beta = simParams->gbis_beta;
00321 BigReal gamma = simParams->gbis_gamma;
00322 BigReal alphaMax = simParams->alpha_max;
00323 int exclude = simParams->exclude;
00324 BigReal scale14 = (exclude==SCALED14) ? simParams->scale14 : 1.0;
00325
00326
00327 #if PRINT > 0
00328
00329 #endif
00330
00331
00332 BigReal kappa = simParams->kappa;
00333
00334
00335 gbInterEnergy = 0;
00336 gbSelfEnergy = 0;
00337 coulEnergy = 0;
00338 BigReal fx, fy, fz;
00339 loop1Iter = 0;
00340 loop2Iter = 0;
00341 loop3Iter = 0;
00342 int loop1flops = 0;
00343 int loop2flops = 0;
00344 int loop3flops = 0;
00345 BigReal loop1time = 0;
00346 BigReal loop2time = 0;
00347 BigReal loop3time = 0;
00348
00349 Position ri, rj;
00350 BigReal rhoi, rhoi0, rhois;
00351 BigReal rhoj, rhoj0, rhojs;
00352 BigReal dx, dy, dz;
00353 BigReal r, r2, r_i, qiqj;
00354 BigReal rnx, rny, rnz;
00355 BigReal bornRad;
00356 BigReal aiaj, aiaj4, expr2aiaj4, fij, finv, expkappa, Dij;
00357 BigReal rc, shift, ddrshift;
00358 BigReal tmp_dEda, dEdai, dEdaj, gbEij;
00359 BigReal ddrfij, ddrfinv, ddrDij, ddrGbEij;
00360 BigReal coulEij, ddrCoulEij;
00361 BigReal forceCoul, forcedEdr, forceAlpha;
00362 BigReal hij, hji, ddrHij, ddrHji;
00363 BigReal tanhi, daidr, dajdr;
00364 BigReal nbetapsi, gammapsi2;
00365 vect *coulForce = new vect[numAtoms];
00366 vect *gbForceR = new vect[numAtoms];
00367 vect *gbForceA = new vect[numAtoms];
00368 BigReal *a = new BigReal[numAtoms];
00369 BigReal *psi = new BigReal[numAtoms];
00370 BigReal *dEda = new BigReal[numAtoms];
00371 BigReal *ddrHijPrefix = new BigReal[numAtoms];
00372
00373
00374 for (int i = 0; i < numAtoms; i++) {
00375 force[i].force.x =0.0;
00376 force[i].force.y =0.0;
00377 force[i].force.z =0.0;
00378 slowForce[i].force.x =0.0;
00379 slowForce[i].force.y =0.0;
00380 slowForce[i].force.z =0.0;
00381 dEda[i] = 0;
00382
00383 psi[i] = 0;
00384 ddrHijPrefix[i] = 0;
00385 coulForce[i].x = 0;
00386 coulForce[i].y = 0;
00387 coulForce[i].z = 0;
00388 gbForceR[i].x = 0;
00389 gbForceR[i].y = 0;
00390 gbForceR[i].z = 0;
00391 gbForceA[i].x = 0;
00392 gbForceA[i].y = 0;
00393 gbForceA[i].z = 0;
00394 }
00395
00396 t2 = clock();
00397 inittime = (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00398
00399
00400
00401
00402
00403
00404
00405
00406 double missTime = 0;
00407 double hitTime = 0;
00408 t1 = clock();
00409 int dij, dji;
00410 BigReal dhij, dhji;
00411
00412
00413 for (int i = 0; i < numAtoms; i++) {
00414 ri = coord[i].position;
00415 for (int j = i+1; j < numAtoms; j++) {
00416 rj = coord[j].position;
00417 dx = (ri.x - rj.x);
00418 dy = (ri.y - rj.y);
00419 dz = (ri.z - rj.z);
00420 r2 = dx*dx+dy*dy+dz*dz;
00421 if (r2 > a_cut_ps2) continue;
00422 r_i = 1.0/sqrt(r2);;
00423 r = 1.0/r_i;;
00424 loop1Iter++;
00425 BigReal rhoi0 = coord[i].rho0;
00426 BigReal rhojs = coord[j].rhos;
00427 BigReal rhoj0 = coord[j].rho0;
00428 BigReal rhois = coord[i].rhos;
00429 CalcHPairSer(r,r2,r_i,a_cut, coord[i].rho0, coord[j].rhos,
00430 coord[j].rho0, coord[i].rhos,dij,dji,hij,hji);
00431
00432 #ifdef PRINT_COMP
00433 CkPrintf("PSI(%04i)[%04i,%04i] = %i%i % .4e % .4e\n",timestep,coord[i].id,coord[j].id,dij,dji,hij,hji);
00434
00435
00436 #endif
00437 psi[i] += hij;
00438 psi[j] += hji;
00439
00440 }
00441 }
00442 t2 = clock();
00443 psitime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00444 loop1time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00445 t1 = clock();
00446
00447
00448
00449
00450
00451
00452
00453
00454 BigReal totPsi = 0;
00455 for (int i = 0; i < numAtoms; i++) {
00456
00457
00458 rhoi0 = coord[i].rho0;
00459 rhoi = coord[i].rho;
00460
00461 totPsi += psi[i];
00462 psi[i] *= rhoi0;
00463 bornRad=1/(1/rhoi0-1/rhoi*tanh(psi[i]*(delta+psi[i]*(-beta+gamma*psi[i]))));
00464 bornRad = (1.0/bornRad < 1.0/alphaMax) ? alphaMax : bornRad;
00465 a[i] = bornRad;
00466 #ifdef PRINT_COMP
00467 CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",timestep,coord[i].id,bornRad);
00468 #endif
00469 }
00470 t2 = clock();
00471 alphatime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00472 t1 = clock();
00473
00474
00475
00476
00477
00478
00479
00480
00481 for (int i = 0; i < numAtoms; i++) {
00482 ri = coord[i].position;
00483 for (int j = i+1; j < numAtoms; j++) {
00484 rj = coord[j].position;
00485 dx = (ri.x - rj.x);
00486 dy = (ri.y - rj.y);
00487 dz = (ri.z - rj.z);
00488 r2 = dx*dx+dy*dy+dz*dz;
00489 if (r2 > cutoff2) continue;
00490
00491
00492 loop2Iter++;
00493 qiqj = (-COUL_CONST*coord[i].charge);
00494 qiqj *= coord[j].charge;
00495 r_i = 1.0/sqrt(r2);
00496 r = 1.0/r_i;
00497 rnx = dx*r_i;
00498 rny = dy*r_i;
00499 rnz = dz*r_i;
00500
00501
00502 BigReal scale, ddrScale;
00503
00504 if ( false ) {
00505 BigReal ratio = r2 / cutoff2;
00506 scale = ratio - 1;
00507 scale *= scale;
00508 ddrScale = 4.0*r*(r2-cutoff2)/cutoff2/cutoff2;
00509 } else {
00510 scale = 1;
00511 ddrScale = 0;
00512 }
00513
00514
00515
00516
00517 int exclij = 0;
00518
00519
00520
00521 Phase2_PairSer(r,r2,r_i,qiqj,a[i],a[j],epsilon_p_i,epsilon_s_i,kappa,exclij,
00522 scale14,stat,coulEij,ddrCoulEij,gbEij,ddrGbEij,dEdai,dEdaj);
00523 int id1 = coord[i].id;
00524 int id2 = coord[j].id;
00525 if (id1 > id2 ) {
00526 int tmp = id2;
00527 id2 = id1;
00528 id1 = tmp;
00529 }
00530
00531
00532 gbInterEnergy += gbEij *scale;
00533 forcedEdr = -(ddrGbEij)*scale-(gbEij)*ddrScale;
00534
00535
00536 fx = rnx*forcedEdr;
00537 fy = rny*forcedEdr;
00538 fz = rnz*forcedEdr;
00539 #ifdef PRINT_COMP
00540 CkPrintf("DEDR(%04i)[%04i,%04i] = % .4e\n",timestep,i,j,forcedEdr);
00541 CkPrintf("DASM(%04i)[%04i,%04i] = % .4e % .4e\n",timestep,i,j,dEdai*scale,dEdaj*scale);
00542 CkPrintf("P2RM(%04i)[%04i,%04i] = % .4e % .4e % .4e % .4e % .4e\n",timestep,i,j, r, a[i],a[j],epsilon_p_i,epsilon_s_i,kappa);
00543
00544 #endif
00545 gbForceR[i].x += fx;
00546 gbForceR[i].y += fy;
00547 gbForceR[i].z += fz;
00548 gbForceR[j].x -= fx;
00549 gbForceR[j].y -= fy;
00550 gbForceR[j].z -= fz;
00551
00552
00553 if (stat) {
00554 dEda[i] += dEdai*scale;
00555 dEda[j] += dEdaj*scale;
00556 }
00557
00558 }
00559 }
00560 t2 = clock();
00561 dEdasumtime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00562 loop2time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00563 t1 = clock();
00564
00565
00566
00567
00568
00569
00570
00571
00572 for (int i = 0; i < numAtoms; i++) {
00573
00574 fij = a[i];
00575 expkappa = exp(-kappa*fij);
00576 Dij = epsilon_p_i - expkappa*epsilon_s_i;
00577
00578 gbEij = -COUL_CONST*coord[i].charge*coord[i].charge*Dij/fij;
00579 gbSelfEnergy += 0.5*gbEij;
00580
00581
00582
00583
00584 if (stat) {
00585 dEdai = -0.5*COUL_CONST*coord[i].charge*coord[i].charge
00586 *(kappa*epsilon_s_i*expkappa-Dij/fij)/a[i];
00587
00588 dEda[i] += dEdai;
00589 BigReal dedasum = dEda[i];
00590
00591
00592 rhoi0 = coord[i].rho0;
00593 rhoi = rhoi0+rho_0;
00594 BigReal psii = psi[i];
00595 nbetapsi = -beta*psii;
00596 gammapsi2 = gamma*psii*psii;
00597 tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
00598 daidr = a[i]*a[i]*rhoi0/rhoi*(1-tanhi*tanhi)
00599 * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
00600 ddrHijPrefix[i] = daidr*dEda[i];
00601 #ifdef PRINT_COMP
00602 CkPrintf("DHDR(%04i)[%04i] = % .4e\n",timestep,coord[i].id, ddrHijPrefix[i]);
00603 #endif
00604 }
00605 }
00606 t2 = clock();
00607 dEdaprefixtime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00608
00609 t1 = clock();
00610
00611
00612
00613
00614
00615
00616 if (stat) {
00617 BigReal dhij, dhji;
00618 int dij, dji;
00619
00620 for (int i = 0; i < numAtoms; i++) {
00621 ri = coord[i].position;
00622
00623
00624 for (int j = i+1; j < numAtoms; j++) {
00625 rj = coord[j].position;
00626 dx = (ri.x - rj.x);
00627 dy = (ri.y - rj.y);
00628 dz = (ri.z - rj.z);
00629 r2 = dx*dx+dy*dy+dz*dz;
00630
00631 if (r2 > a_cut_ps2) continue;
00632 r_i = 1.0/sqrt(r2);
00633 loop3Iter++;
00634
00635 r = 1.0/r_i;
00636 CalcDHPairSer(r,r2,r_i,a_cut,
00637 coord[i].rho0,coord[j].rhos,
00638 coord[j].rho0,coord[i].rhos,
00639 dij,dji,dhij,dhji);
00640
00641
00642 forceAlpha = -r_i*(ddrHijPrefix[i]*dhij+ddrHijPrefix[j]*dhji);
00643 #ifdef PRINT_COMP
00644 CkPrintf("DEDA(%04i)[%04i,%04i] = %i%i % .4e % .4e % .4e\n",timestep,i,j,dij,dji,dhij,dhji, forceAlpha/r_i);
00645 #endif
00646 fx = dx * forceAlpha;
00647 fy = dy * forceAlpha;
00648 fz = dz * forceAlpha;
00649
00650 gbForceA[i].x += fx;
00651 gbForceA[i].y += fy;
00652 gbForceA[i].z += fz;
00653 gbForceA[j].x -= fx;
00654 gbForceA[j].y -= fy;
00655 gbForceA[j].z -= fz;
00656 }
00657 }
00658 }
00659 t2 = clock();
00660 dEdalooptime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00661 loop3time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00662 t1 = clock();
00663
00664
00665 for (int i = 0; i < numAtoms; i++) {
00666 force[i].force.x =0.0;
00667 force[i].force.y =0.0;
00668 force[i].force.z =0.0;
00669 slowForce[i].force.x =0.0;
00670 slowForce[i].force.y =0.0;
00671 slowForce[i].force.z =0.0;
00672
00673 #if GBIS_COUL_FORCE
00674 force[i].force.x += coulForce[i].x;
00675 force[i].force.y += coulForce[i].y;
00676 force[i].force.z += coulForce[i].z;
00677 #endif
00678
00679 #if GBIS_DEDR_FORCE
00680 force[i].force.x += gbForceR[i].x;
00681 force[i].force.y += gbForceR[i].y;
00682 force[i].force.z += gbForceR[i].z;
00683 #endif
00684
00685 #if GBIS_DEDA_FORCE
00686 if (stat) {
00687 slowForce[i].force.x += gbForceA[i].x;
00688 slowForce[i].force.y += gbForceA[i].y;
00689 slowForce[i].force.z += gbForceA[i].z;
00690
00691 }
00692
00693
00694
00695 #endif
00696
00697
00698 #ifdef PRINT_SERFORCES
00699 BigReal fC, fR, fA;
00700 fx = 0;
00701 fy = 0;
00702 fz = 0;
00703 #if GBIS_COUL_FORCE
00704 fx += coulForce[i].x;
00705 fy += coulForce[i].y;
00706 fz += coulForce[i].z;
00707
00708
00709 #endif
00710 #if GBIS_DEDR_FORCE
00711 fx += gbForceR[i].x;
00712 fy += gbForceR[i].y;
00713 fz += gbForceR[i].z;
00714
00715
00716 #endif
00717 #if GBIS_DEDA_FORCE
00718 fx += gbForceA[i].x;
00719 fy += gbForceA[i].y;
00720 fz += gbForceA[i].z;
00721
00722
00723 #endif
00724
00725 #endif //if print>1
00726 }
00727 t2 = clock();
00728 inittime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00729
00730 timestep++;
00731
00732 delete[] coulForce;
00733 delete[] gbForceR;
00734 delete[] gbForceA;
00735 delete[] a;
00736 delete[] psi;
00737 delete[] dEda;
00738 delete[] ddrHijPrefix;
00739
00740 }
00741
00742 void ComputeGBISserMgr::recvCoord(GBISCoordMsg *msg) {
00743
00744 if ( ! numSources ) {
00745 numSources = (PatchMap::Object())->numNodesWithPatches();
00746 coordMsgs = new GBISCoordMsg*[numSources];
00747 for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
00748 numArrived = 0;
00749 numAtoms = Node::Object()->molecule->numAtoms;
00750 coord = new ComputeGBISAtom[numAtoms];
00751 force = new ExtForce[numAtoms];
00752 slowForce = new ExtForce[numAtoms];
00753 timestep = msg[0].sequence;
00754 }
00755
00756
00757
00758 int i;
00759 for ( i=0; i < msg->numAtoms; ++i ) {
00760 coord[msg->coord[i].id] = msg->coord[i];
00761 }
00762
00763 coordMsgs[numArrived] = msg;
00764 ++numArrived;
00765
00766 if ( numArrived < numSources ) return;
00767 numArrived = 0;
00768
00769
00770
00771
00772
00773
00774
00775 t_start = clock();
00776
00777 double namdtime = (double)((double)t_start-(double)t_stop)/CLOCKS_PER_SEC;
00778 totalnamdtime += namdtime;
00779
00780 calcGBISReg(msg->doSlow);
00781
00782 t_stop = clock();
00783 double gbistime = (double)((double)t_stop-(double)t_start)/CLOCKS_PER_SEC;
00784
00785
00786 totalgbistime += gbistime;
00787
00788 if (timestep % BENCH_PERIOD == 0) {
00789 printf("\n");
00790 printf("GBIS: t_GB=%f sec for %i steps\n",totalgbistime,BENCH_PERIOD);
00791 printf("GBIS: t_MD=%f sec for %i steps\n",totalnamdtime,BENCH_PERIOD);
00792 printf("GBIS: init=%f sec for %i steps\n", inittime, BENCH_PERIOD);
00793 printf("GBIS: psi=%f sec for %i steps\n", psitime, BENCH_PERIOD);
00794 printf("GBIS: alpha=%f sec for %i steps\n", alphatime, BENCH_PERIOD);
00795 printf("GBIS: dEdasum=%f sec for %i steps\n", dEdasumtime, BENCH_PERIOD);
00796 printf("GBIS: dEdaprefix=%f sec for %i steps\n",dEdaprefixtime,BENCH_PERIOD);
00797 printf("GBIS: dEdaloop=%f sec for %i steps\n", dEdalooptime,BENCH_PERIOD);
00798 printf("GBIS: loop1=%i iters\n", loop1Iter);
00799 printf("GBIS: loop2=%i iters\n", loop2Iter);
00800 printf("GBIS: loop3=%i iters\n", loop3Iter);
00801 printf("\n");
00802 totalgbistime = 0;
00803 totalnamdtime = 0;
00804 inittime = 0;
00805 psitime = 0;
00806 alphatime = 0;
00807 dEdrtime = 0;
00808 dEdasumtime = 0;
00809 dEdaprefixtime = 0;
00810 dEdalooptime = 0;
00811 }
00812
00813
00814
00815 for ( int j=0; j < numSources; ++j ) {
00816 GBISCoordMsg *cmsg = coordMsgs[j];
00817 coordMsgs[j] = 0;
00818 GBISForceMsg *fmsg = new (cmsg->numAtoms, cmsg->numAtoms, 0) GBISForceMsg;
00819
00820 for ( int i=0; i < cmsg->numAtoms; ++i ) {
00821 fmsg->force[i] = force[cmsg->coord[i].id];
00822 fmsg->slowForce[i] = slowForce[cmsg->coord[i].id];
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835 }
00836
00837
00838 if ( ! j ) {
00839 #if GBIS_DEDR_FORCE
00840 fmsg->gbSelfEnergy = gbSelfEnergy;
00841 fmsg->gbInterEnergy = gbInterEnergy;
00842 #else
00843 fmsg->gbSelfEnergy = 0;
00844 fmsg->gbInterEnergy = 0;
00845 #endif
00846 #if GBIS_COUL_FORCE
00847 fmsg->coulEnergy = coulEnergy;
00848 #else
00849 fmsg->coulEnergy = 0;
00850 #endif
00851 } else {
00852 fmsg->gbSelfEnergy = 0;
00853 fmsg->gbInterEnergy = 0;
00854 fmsg->coulEnergy = 0;
00855 }
00856 gbisProxy[cmsg->sourceNode].recvForce(fmsg);
00857 delete cmsg;
00858 }
00859
00860 }
00861
00862 void ComputeGBISserMgr::recvForce(GBISForceMsg *msg) {
00863
00864 gbisCompute->saveResults(msg);
00865 delete oldmsg;
00866 oldmsg = msg;
00867 }
00868
00869 void ComputeGBISser::saveResults(GBISForceMsg *msg) {
00870
00871 ResizeArrayIter<PatchElem> ap(patchList);
00872
00873 ExtForce *results_ptr = msg->force;
00874 ExtForce *results_ptr_slow = msg->slowForce;
00875
00876
00877 for (ap = ap.begin(); ap != ap.end(); ap++) {
00878
00879 Results *r = (*ap).forceBox->open();
00880 Force *f = r->f[Results::nbond];
00881 Force *sf = r->f[Results::slow];
00882 int numAtoms = (*ap).p->getNumAtoms();
00883
00884 for(int i=0; i<numAtoms; ++i) {
00885 f[i] += results_ptr->force;
00886 sf[i] += results_ptr_slow->force;
00887
00888 ++results_ptr;
00889 ++results_ptr_slow;
00890 }
00891
00892 (*ap).forceBox->close(&r);
00893 }
00894
00895
00896 reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbInterEnergy;
00897 reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbSelfEnergy;
00898
00899 reduction->submit();
00900 }
00901
00902
00903
00904
00905 void ComputeGBISser::doWork()
00906 {
00907
00908 ResizeArrayIter<PatchElem> ap(patchList);
00909
00910 #if 1
00911
00912 if ( ! patchList[0].p->flags.doNonbonded )
00913 {
00914 for (ap = ap.begin(); ap != ap.end(); ap++) {
00915 CompAtom *x = (*ap).positionBox->open();
00916
00917 Results *r = (*ap).forceBox->open();
00918 (*ap).positionBox->close(&x);
00919
00920 (*ap).forceBox->close(&r);
00921 }
00922 reduction->submit();
00923 return;
00924 }
00925 #endif
00926
00927
00928
00929 int numLocalAtoms = 0;
00930 for (ap = ap.begin(); ap != ap.end(); ap++) {
00931 numLocalAtoms += (*ap).p->getNumAtoms();
00932 }
00933
00934 GBISCoordMsg *msg = new (numLocalAtoms, 0) GBISCoordMsg;
00935 msg->sourceNode = CkMyPe();
00936 msg->numAtoms = numLocalAtoms;
00937 ComputeGBISAtom *data_ptr = msg->coord;
00938 SimParameters *simParams = Node::Object()->simParameters;
00939 msg->doSlow = patchList[0].p->flags.doFullElectrostatics;
00940
00941 msg->sequence = sequence();
00942
00943
00944 for (ap = ap.begin(); ap != ap.end(); ap++) {
00945 FullAtomList &atoms = (*ap).p->getAtomList();
00946 CompAtom *x = (*ap).positionBox->open();
00947 CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00948 int numAtoms = (*ap).p->getNumAtoms();
00949 for(int i=0; i<numAtoms; ++i)
00950 {
00951 data_ptr->position = x[i].position;
00952 data_ptr->charge = x[i].charge;
00953 data_ptr->mass = atoms[i].mass;
00954 data_ptr->id = xExt[i].id;
00955 data_ptr->rho = MassToRadius(data_ptr->mass);
00956 SimParameters *simParams = Node::Object()->simParameters;
00957 data_ptr->rho0 = data_ptr->rho - simParams->coulomb_radius_offset;
00958 data_ptr->rhos = data_ptr->rho0 * MassToScreen(data_ptr->mass);
00959 data_ptr->vdwType = x[i].vdwType;
00960 ++data_ptr;
00961 }
00962
00963 #if 0
00964 if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00965 else { (*ap).positionBox->close(&x); }
00966 #endif
00967 (*ap).positionBox->close(&x);
00968 }
00969
00970 CProxy_ComputeGBISserMgr gbisProxy(CkpvAccess(BOCclass_group).computeGBISserMgr);
00971 gbisProxy[0].recvCoord(msg);
00972
00973 }
00974 #include "ComputeGBISserMgr.def.h"
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987 inline void h0Ser ( BigReal r, BigReal r2, BigReal ri,
00988 Real rc, BigReal r0, BigReal rs, BigReal & h ) {
00989 h = 0;
00990 }
00991 inline void dh0Ser ( BigReal r, BigReal r2, BigReal ri,
00992 Real rc, BigReal r0, BigReal rs, BigReal & dh ) {
00993 dh = 0;
00994 }
00995
00996 inline void h1Ser ( BigReal r, BigReal r2, BigReal ri,
00997 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
00998
00999 BigReal rci = 1.0/rc;
01000 BigReal rmrs = r-rs;
01001 BigReal rmrsi = 1.0/rmrs;
01002 BigReal rmrs2 = rmrs*rmrs;
01003 BigReal rs2 = rs*rs;
01004 BigReal logr = log(rmrs*rci);
01005 BigReal rci2 = rci*rci;
01006 h = 0.125*ri*(1 + 2*r*rmrsi + rci2*(r2 - 4*rc*r - rs2) + 2*logr);
01007 }
01008 inline void dh1Ser ( BigReal r, BigReal r2, BigReal ri,
01009 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01010
01011 BigReal rci = 1.0/rc;
01012 BigReal rmrs = r-rs;
01013 BigReal rmrsi = 1.0/rmrs;
01014 BigReal rmrs2 = rmrs*rmrs;
01015 BigReal rs2 = rs*rs;
01016 BigReal logr = log(rmrs*rci);
01017 BigReal rci2 = rci*rci;
01018 dh = ri*ri*(-0.25*logr - (rc*rc - rmrs2)*(rs2 + r2)*0.125*rci2*rmrsi*rmrsi);
01019 }
01020
01021 inline void h2Ser ( BigReal r, BigReal r2, BigReal ri,
01022 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01023
01024 BigReal k = rs*ri; k*=k;
01025 h = rs*ri*ri*k*(TA+k*(TB+k*(TC+k*(TD+k*TE))));
01026 }
01027 inline void dh2Ser ( BigReal r, BigReal r2, BigReal ri,
01028 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01029
01030 BigReal k = rs*ri; k*=k;
01031 dh = -rs*ri*ri*ri*k*(DA+k*(DB+k*(DC+k*(DD+k*DE))));
01032 }
01033
01034 inline void h3Ser ( BigReal r, BigReal r2, BigReal ri,
01035 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01036 BigReal r2mrs2i = 1.0/(r2-rs*rs);
01037 h = 0.5 * ( rs*r2mrs2i + 0.5 * log((r-rs)/(r+rs))*ri );
01038 }
01039 inline void dh3Ser ( BigReal r, BigReal r2, BigReal ri,
01040 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01041 BigReal rs2 = rs*rs;
01042 BigReal r2mrs2i = 1.0/(r2-rs2);
01043 dh = -0.25*ri*(2*(r2+rs2)*rs*r2mrs2i*r2mrs2i + ri*log((r-rs)/(r+rs)));
01044 }
01045
01046 inline void h4Ser ( BigReal r, BigReal r2, BigReal ri,
01047 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01048 BigReal ri2 = ri*ri;
01049 BigReal r02 = r0*r0;
01050 BigReal rs2 = rs*rs;
01051 BigReal r0i = 1.0/r0;
01052 BigReal rspri = 1.0/(r+rs);
01053 BigReal logr = log(r0*rspri);
01054 BigReal r02mrs2 = r02-rs2;
01055 BigReal rilogr = ri*logr;
01056 h = 0.25*( r0i*(2-0.5*(r0i*ri*(r2 + r02 - rs2))) - rspri + rilogr );
01057 }
01058 inline void dh4Ser ( BigReal r, BigReal r2, BigReal ri,
01059 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01060 BigReal ri2 = ri*ri;
01061 BigReal r02 = r0*r0;
01062 BigReal rs2 = rs*rs;
01063 BigReal r0i = 1.0/r0;
01064 BigReal rspri = 1.0/(r+rs);
01065 BigReal logr = log(r0*rspri);
01066 BigReal r02mrs2 = r02-rs2;
01067 BigReal rilogr = ri*logr;
01068 dh = 0.25*( (-0.5+(r2*r02mrs2 - 2*r*rs*rs2+rs2*r02mrs2)
01069 * 0.5*ri2*rspri*rspri)*r0i*r0i - ri*rilogr );
01070 }
01071
01072 inline void h5Ser ( BigReal r, BigReal r2, BigReal ri,
01073 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01074 BigReal rs2 = rs*rs;
01075 BigReal r2mrs2i = 1/(r2-rs2);
01076 BigReal rsr2mrs2i = rs*r2mrs2i;
01077 BigReal rprs = r+rs;
01078 BigReal rmrs = r-rs;
01079 BigReal logr = 0.5*ri*log(-rmrs/rprs);
01080 h = 0.5*( rsr2mrs2i + 2/r0 + logr );
01081 }
01082 inline void dh5Ser ( BigReal r, BigReal r2, BigReal ri,
01083 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01084 BigReal rs2 = rs*rs;
01085 BigReal r2mrs2i = 1/(r2-rs2);
01086 BigReal rsr2mrs2i = rs*r2mrs2i;
01087 BigReal rprs = r+rs;
01088 BigReal rmrs = r-rs;
01089 BigReal logr = 0.5*ri*log(-rmrs/rprs);
01090 dh = -0.5*ri*((rs2+r2)*rsr2mrs2i*r2mrs2i+logr );
01091 }
01092
01093 inline void h6Ser ( BigReal r, BigReal r2, BigReal ri,
01094 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01095 h = 0;
01096 }
01097 inline void dh6Ser ( BigReal r, BigReal r2, BigReal ri,
01098 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01099 dh = 0;
01100 }
01101
01102 inline void CalcHSer ( BigReal r, BigReal r2, BigReal ri,
01103 BigReal rc, BigReal r0, BigReal rs, BigReal & h, int & d) {
01104 if( r <= rc - rs && r > 4*rs ) {
01105 h2Ser(r,r2,ri,rc,r0,rs,h); d = 2;
01106 } else if (r <= rc + rs && r > rc - rs) {
01107 h1Ser(r,r2,ri,rc,r0,rs,h); d = 1;
01108 } else if (r > rc + rs) {
01109 h0Ser(r,r2,ri,rc,r0,rs,h); d = 0;
01110 } else if( r <= 4*rs && r > r0 + rs ) {
01111 h3Ser(r,r2,ri,rc,r0,rs,h); d = 3;
01112 } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {
01113 h4Ser(r,r2,ri,rc,r0,rs,h); d = 4;
01114 } else if (r0 < rs ) {
01115 h5Ser(r,r2,ri,rc,r0,rs,h); d = 5;
01116 } else {
01117 h6Ser(r,r2,ri,rc,r0,rs,h); d = 6;
01118 }
01119 }
01120 inline void CalcDHSer ( BigReal r, BigReal r2, BigReal ri,
01121 BigReal rc, BigReal r0, BigReal rs, BigReal & dh, int & d) {
01122 if( r <= rc - rs && r > 4*rs ) {
01123 dh2Ser(r,r2,ri,rc,r0,rs,dh); d = 2;
01124 } else if (r <= rc + rs && r > rc - rs) {
01125 dh1Ser(r,r2,ri,rc,r0,rs,dh); d = 1;
01126 } else if (r > rc + rs) {
01127 dh0Ser(r,r2,ri,rc,r0,rs,dh); d = 0;
01128 } else if( r <= 4*rs && r > r0 + rs ) {
01129 dh3Ser(r,r2,ri,rc,r0,rs,dh); d = 3;
01130 } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {
01131 dh4Ser(r,r2,ri,rc,r0,rs,dh); d = 4;
01132 } else if (r0 < rs ) {
01133 dh5Ser(r,r2,ri,rc,r0,rs,dh); d = 5;
01134 } else {
01135 dh6Ser(r,r2,ri,rc,r0,rs,dh); d = 6;
01136 }
01137 }
01138 inline void CalcHPairSer (
01139 BigReal r,
01140 BigReal r2,
01141 BigReal ri,
01142 BigReal rc,
01143 BigReal ri0,
01144 BigReal rjs,
01145 BigReal rj0,
01146 BigReal ris,
01147 int & dij,
01148 int & dji,
01149 BigReal & hij,
01150 BigReal & hji
01151 ) {
01152 CalcHSer(r,r2,ri,rc,ri0,rjs,hij,dij);
01153 CalcHSer(r,r2,ri,rc,rj0,ris,hji,dji);
01154 }
01155 inline void CalcDHPairSer (
01156 BigReal r,
01157 BigReal r2,
01158 BigReal ri,
01159 BigReal rc,
01160 BigReal ri0,
01161 BigReal rjs,
01162 BigReal rj0,
01163 BigReal ris,
01164 int & dij,
01165 int & dji,
01166 BigReal & dhij,
01167 BigReal & dhji
01168 ) {
01169 CalcDHSer(r,r2,ri,rc,ri0,rjs,dhij,dij);
01170 CalcDHSer(r,r2,ri,rc,rj0,ris,dhji,dji);
01171 }
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183 inline void Calc_Coul_PairSer(
01184 BigReal r_i,
01185 BigReal qiqj,
01186 BigReal epsilon_p_i,
01187 int exclij,
01188 BigReal scale14,
01189 BigReal & coulE,
01190 BigReal & ddrCoulE
01191 ) {
01192 if (exclij != EXCHCK_FULL) {
01193
01194 coulE = -qiqj*epsilon_p_i*r_i;
01195
01196
01197 if (exclij == EXCHCK_MOD)
01198 coulE *= scale14;
01199 ddrCoulE = -r_i*coulE;
01200 } else {
01201 coulE = 0;
01202 ddrCoulE = 0;
01203 }
01204 }
01205
01206
01207
01208
01209
01210 inline void Calc_dEdr_PairSer(
01211 BigReal r,
01212 BigReal r2,
01213 BigReal qiqj,
01214 BigReal ai,
01215 BigReal aj,
01216 BigReal kappa,
01217 BigReal epsilon_p_i,
01218 BigReal epsilon_s_i,
01219 BigReal & aiaj,
01220 BigReal & expr2aiaj4,
01221 BigReal & fij,
01222 BigReal & f_i,
01223 BigReal & expkappa,
01224 BigReal & Dij,
01225 BigReal & gbE,
01226 BigReal & ddrGbE
01227 ) {
01228
01229 BigReal aiaj4,ddrDij,ddrf_i,ddrfij;
01230
01231
01232 aiaj = ai*aj;
01233 aiaj4 = 4*aiaj;
01234 expr2aiaj4 = exp(-r2/aiaj4);
01235 fij = sqrt(r2+aiaj*expr2aiaj4);
01236 f_i = 1/fij;
01237 if (kappa > 0)
01238 expkappa = exp(-kappa*fij);
01239 else
01240 expkappa = 1.0;
01241 Dij = epsilon_p_i - expkappa*epsilon_s_i;
01242 gbE = qiqj*Dij*f_i;
01243
01244
01245 ddrfij = r*f_i*(1 - 0.25*expr2aiaj4);
01246 ddrf_i = -ddrfij*f_i*f_i;
01247 ddrDij = kappa*expkappa*ddrfij*epsilon_s_i;
01248 ddrGbE = qiqj*(ddrDij*f_i+Dij*ddrf_i);
01249 }
01250
01251
01252
01253
01254
01255 inline void Calc_dEda_PairSer(
01256 BigReal r2,
01257 BigReal ai,
01258 BigReal aj,
01259 BigReal qiqj,
01260 BigReal kappa,
01261 BigReal aiaj,
01262 BigReal expkappa,
01263 BigReal expr2aiaj4,
01264 BigReal fij,
01265 BigReal f_i,
01266 BigReal Dij,
01267 BigReal epsilon_s_i,
01268 BigReal & dEdai,
01269 BigReal & dEdaj
01270 ) {
01271
01272 BigReal tmp_dEda = 0.5*qiqj*f_i*f_i
01273 *(kappa*epsilon_s_i*expkappa-Dij*f_i)
01274 *(aiaj+0.25*r2)*expr2aiaj4;
01275 dEdai = tmp_dEda/ai;
01276 dEdaj = tmp_dEda/aj;
01277 }
01278
01279
01280
01281
01282
01283 inline void Phase2_PairSer(
01284
01285
01286 BigReal r,
01287 BigReal r2,
01288 BigReal r_i,
01289 BigReal qiqj,
01290 BigReal ai,
01291 BigReal aj,
01292 BigReal epsilon_p_i,
01293 BigReal epsilon_s_i,
01294 BigReal kappa,
01295 int exclij,
01296 BigReal scale14,
01297 int doSlow,
01298
01299
01300 BigReal & coulEij,
01301 BigReal & ddrCoulEij,
01302 BigReal & gbEij,
01303 BigReal & ddrGbEij,
01304 BigReal & dEdai,
01305 BigReal & dEdaj
01306 ) {
01307
01308
01309
01310 coulEij = 0;
01311 ddrCoulEij = 0;
01312
01313
01314 BigReal aiaj,expr2aiaj4,fij,f_i,expkappa,Dij;
01315 Calc_dEdr_PairSer(r,r2,qiqj,ai,aj,kappa,epsilon_p_i,epsilon_s_i,
01316 aiaj,expr2aiaj4,fij,f_i,expkappa,Dij,gbEij,ddrGbEij);
01317
01318
01319 if (doSlow) {
01320 Calc_dEda_PairSer(r2,ai,aj,qiqj,kappa,aiaj,expkappa,expr2aiaj4,
01321 fij,f_i,Dij,epsilon_s_i,dEdai,dEdaj);
01322 }
01323 }
01324