NAMD: ComputeGBISser.C Source File

00001 
00007 /*******************************************************************************
00008  *******************************************************************************
00009   This serial version of GBIS is out of date and hasn't been vetted;
00010   it is to be used as a way of testing new implicit solvent models
00011   since all atomic coordinates are gathered into a single Compute.
00012  *******************************************************************************
00013  ******************************************************************************/
00014 
00015 /*
00016   print debug
00017   1 - once per step
00018   2 - once per atom
00019   3 - once per pair
00020 */
00021 //#define PRINT_SERFORCES
00022 #define BENCH_PERIOD 1000
00023 
00024 //sum all forces
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 // #define DEBUGM
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(//doesn't do self energies
00076 
00077 //input values
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 //return values
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,//distance
00102   BigReal r2,
00103   BigReal ri,
00104   BigReal rc,//cutoff
00105   BigReal ri0,
00106   BigReal rjs,
00107   BigReal rj0,
00108   BigReal ris,
00109   int & dij,//domain 1
00110   int & dji,//domain 2
00111   BigReal & dhij,
00112   BigReal & dhji);
00113 inline void CalcDHPairSer (
00114   BigReal r,//distance
00115   BigReal r2,
00116   BigReal ri,
00117   BigReal rc,//cutoff
00118   BigReal ri0,
00119   BigReal rjs,
00120   BigReal rj0,
00121   BigReal ris,
00122   int & dij,//domain 1
00123   int & dji,//domain 2
00124   BigReal & dhij,
00125   BigReal & dhji);
00126 
00127 
00128 struct ComputeGBISAtom {
00129     Position position;
00130     Real charge;
00131     Real rho;//coulomb radius
00132     Real rho0;//coulomb radius
00133     Real rhos;//coulomb radius
00134     Mass mass;
00135     //int type;//atom type for S scaling table
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  * scale from 1 to 0 between transition and cutoff
00248  * s = 1 for [0,trans]
00249  * s = s(r) for [trans, cutoff]
00250  * s = 0 for [cutoff, infinity]
00251  * meets the following 4 conditions:
00252  * (1) s(trans)   = 1
00253  * (3) s'(trans)  = 0
00254  * (2) s(cutoff)  = 0
00255  * (4) s'(cutoff) = 0
00256  ******************************************************************************/
00257 inline void CalcScaleSer (
00258   BigReal r,//rij
00259   BigReal t,//transition distance
00260   BigReal c,//cutoff distance
00261   BigReal & s,//output s(r)
00262   BigReal & d ) {//output s'(r)
00263 
00264   if (r <= t) { //[0,trans]
00265     s = 1;
00266     d = 0;
00267   } else if (r < c) { //[trans,cutoff]
00268     //precompute
00269     BigReal ct = (c-t);
00270     BigReal ct_i = 1/ct;
00271     BigReal ct_i2 = ct_i*ct_i;
00272     //BigReal ct_i4 = ct_i2*ct_i2;
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 { //[cutoff,infinity]
00283     s = 0;
00284     d = 0;
00285   }
00286 }
00287 
00288 
00289 /**********************************************************
00290  *
00291  *  Regular 3 sets of nested loops
00292  *
00293  *  supports 1-2-4 scheme
00294  *
00295  **********************************************************/
00296 void ComputeGBISserMgr::calcGBISReg(int stat) {
00297   t1 = clock();
00298   //printf("GBIS: ComputeGBIS::serial\n");
00299 
00300   /*********************************************************
00301    *  Initializations
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;// + 1.8*0.96;//+max screen radius
00311   BigReal trans = simParams->switchingDist;
00312   BigReal fsMax = 1.728;
00313   BigReal a_cut = simParams->alpha_cutoff-fsMax;//+max screen radius; partial sphere
00314   //a_cut = 10;
00315   BigReal a_cut2 = a_cut*a_cut;
00316   BigReal a_cut_ps = a_cut + fsMax;//+max screen radius; partial sphere
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   //printf("cut=%e, rgbmax=%e, fsmax=%e, gbalpha=%e, gbbeta=%e, gbgamma=%e\n", cutoff, alphaMax, fsMax, delta, beta, gamma);
00329 #endif
00330 
00331   //calc kappa
00332   BigReal kappa = simParams->kappa;
00333 
00334   //debugging variables
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];//dEdr
00367   vect *gbForceA = new vect[numAtoms];//dEda*dadr
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   //init arrays
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;//dEdr
00389     gbForceR[i].y = 0;//dEdr
00390     gbForceR[i].z = 0;//dEdr
00391     gbForceA[i].x = 0;//dEda*dadr
00392     gbForceA[i].y = 0;//dEda*dadr
00393     gbForceA[i].z = 0;//dEda*dadr
00394   }
00395 
00396   t2 = clock();
00397   inittime = (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
00398 
00399   /**********************************************************
00400    *
00401    *  Loop 1 : accumulate Hij into psi array
00402    *  cutoff = a_cut = 25A
00403    *  every 2 timesteps
00404    *
00405    **********************************************************/
00406   double missTime = 0;
00407   double hitTime = 0;
00408   t1 = clock();
00409   int dij, dji;//which hij domain 0123456 - for debugging only
00410   BigReal dhij, dhji;
00411   //int TESTITER = 1;
00412   //for (int k = 0; k < TESTITER; k++)
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       //CkPrintf("S_Hij_%i\n",dij);
00435       //CkPrintf("S_Hji_%i\n",dji);
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    *  Atom 1: Calculate alpha based on phi for each atom
00452    *
00453    **********************************************************/
00454   BigReal totPsi = 0;
00455   for (int i = 0; i < numAtoms; i++) {
00456     //CkPrintf("Srho[%i] = %.3e %.3e\n",coord[i].id,coord[i].rho0,coord[i].rhos);
00457 
00458     rhoi0 = coord[i].rho0;
00459     rhoi = coord[i].rho;
00460     //CkPrintf("GBIS_SER: psi[%03i] = %.4e\n",i,psi[i]);
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    *  Loop 2 : dEda
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);//rptI
00486       dy = (ri.y - rj.y);//rptI
00487       dz = (ri.z - rj.z);//rptI
00488       r2 = dx*dx+dy*dy+dz*dz;//rptI
00489       if (r2 > cutoff2) continue;
00490 
00491       //calculate distance
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;//rptI
00497       rnx = dx*r_i;//normalized vector
00498       rny = dy*r_i;
00499       rnz = dz*r_i;
00500 
00501       //calculate scaling for energy smoothing
00502       BigReal scale, ddrScale;
00503       //CalcScaleSer(r, trans, cutoff, scale, ddrScale);
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       //BigReal ai = a[i];
00516       //BigReal aj = a[j];
00517       int exclij = 0;//molecule->checkexcl(i,j);
00518       //CkPrintf("GBIS_SER_excl[%i,%i] %i\n",i,j,exclij);
00519       //CkPrintf("GBIS_DOFULL = %i\n",stat);
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       }//print ids as lower index, higher index
00530 
00531       //accumulate energies
00532       gbInterEnergy   += gbEij  *scale;
00533       forcedEdr = -(ddrGbEij)*scale-(gbEij)*ddrScale;
00534 
00535       //add GB force
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       //add dEda
00553       if (stat) {
00554         dEda[i] += dEdai*scale;
00555         dEda[j] += dEdaj*scale;
00556       }
00557 
00558     }//end j loop
00559   }//end i loop
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 *   Atom 2 : Calculate dHij Prefix
00570 *
00571 *******************************************************************************/
00572   for (int i = 0; i < numAtoms; i++) {
00573     //add diagonal dEda term
00574     fij = a[i];//inf
00575     expkappa = exp(-kappa*fij);//0
00576     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
00577     //add diagonal GB energy term
00578     gbEij = -COUL_CONST*coord[i].charge*coord[i].charge*Dij/fij;
00579     gbSelfEnergy += 0.5*gbEij;//self energy
00580 //CkPrintf("self_energy[%03i] = %e\n",i,0.5*gbEij);
00581     //CkPrintf("gbSelfEnergy[%03i] = % e\n", coord[i].id,0.5*gbEij);
00582 
00583   //calculate dHij prefix
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       //CkPrintf("SER_dEdai[%03i]%i = % e\n",i,timestep,dEdai);
00588       dEda[i] += dEdai;
00589       BigReal dedasum = dEda[i];
00590       //CkPrintf("SER_dEdaSum[%03i]%i = % e\n",i,timestep,dEda[i]);
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   //gbEnergy = gbInterEnergy + gbSelfEnergy;
00609   t1 = clock();
00610 
00611   /**********************************************************
00612    *
00613    *  Loop 3 : Calculate dEda*dadr Forces
00614    *
00615    **********************************************************/
00616   if (stat) {
00617   BigReal dhij, dhji;
00618   int dij, dji;
00619   //for (int k = 0; k < TESTITER; k++)
00620   for (int i = 0; i < numAtoms; i++) {
00621     ri = coord[i].position;
00622     //CkPrintf("SER_dHdrPrefix[%i]%i = % .5e\n",coord[i].id,timestep,ddrHijPrefix[i]);
00623 
00624     for (int j = i+1; j < numAtoms; j++) {
00625       rj = coord[j].position;
00626       dx = (ri.x - rj.x);//rptI
00627       dy = (ri.y - rj.y);//rptI
00628       dz = (ri.z - rj.z);//rptI
00629       r2 = dx*dx+dy*dy+dz*dz;//rptI
00630       //loop3flops += 9;
00631       if (r2 > a_cut_ps2) continue;
00632       r_i = 1.0/sqrt(r2);//rptI
00633       loop3Iter++;
00634       //calculate dHij for pair
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       //calculate and add dEijdai,j*dai,jdrij force
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     }//end inner summation loop
00657   }//end outer summation loop
00658   }//end if stat
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   //add forces
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       //CkPrintf("SERIAL SLOW %e %e %e\n",gbForceA[i].x,gbForceA[i].y,gbForceA[i].z);
00691     }
00692     //force[i].force.x += gbForceA[i].x;
00693     //force[i].force.y += gbForceA[i].y;
00694     //force[i].force.z += gbForceA[i].z;
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       //fC = sqrt(fx*fx+fy*fy+fz*fz);
00708       //fprintf(stderr, "%i %e %e %e ",i,fx, fy, fz);
00709 #endif
00710 #if GBIS_DEDR_FORCE
00711       fx += gbForceR[i].x;
00712       fy += gbForceR[i].y;
00713       fz += gbForceR[i].z;
00714       //fR = sqrt(fx*fx+fy*fy+fz*fz);
00715       //fprintf(stderr, "%e %e %e",fx, fy, fz);
00716 #endif
00717 #if GBIS_DEDA_FORCE
00718       fx += gbForceA[i].x;
00719       fy += gbForceA[i].y;
00720       fz += gbForceA[i].z;
00721       //fA = sqrt(fx*fx+fy*fy+fz*fz);
00722       //fprintf(stderr, "%5i % .5e % .5e % .5e\n",i, fx, fy, fz);
00723 #endif
00724       //fprintf(stderr, "%5i % .5e % .5e % .5e\n",i, fx, fy, fz);
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   //printf("GBIS calcReg COMPLETE!\n");
00740 }
00741 
00742 void ComputeGBISserMgr::recvCoord(GBISCoordMsg *msg) {
00743   //printf("GBIS recvCoord()\n");
00744   if ( ! numSources ) {//init receiving coord
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   //receive coord
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    *  All sources arrived; calculate energy, forces
00773    *
00774    **********************************************************/
00775   t_start = clock();
00776   //how long did namd take to cycle back here
00777   double namdtime = (double)((double)t_start-(double)t_stop)/CLOCKS_PER_SEC;
00778   totalnamdtime += namdtime;
00779   //choice choose
00780   calcGBISReg(msg->doSlow);
00781   
00782   t_stop = clock();
00783   double gbistime = (double)((double)t_stop-(double)t_start)/CLOCKS_PER_SEC;
00784   //printf("GBIS: elapsednamd(%i)=%f\n", timestep, namdtime);
00785   //printf("GBIS: elapsedgbis(%i)=%f\n", timestep, gbistime);
00786   totalgbistime += gbistime;
00787   //printf("GBIS: total(%i)=%f\n", timestep, totaltime);
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   // distribute forces
00814   //printf("GBIS distributing forces\n");
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     //does this init slowForce, force?
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       if (i == 1000 ) {
00825         printf("%i force: <", i);
00826         printf("%f, ", fmsg->force[i].force.x);
00827         printf("%f, ", fmsg->force[i].force.y);
00828         printf("%f>\n", fmsg->force[i].force.z);
00829         printf("%i slwfr: <", i);
00830         printf("%f, ", fmsg->slowForce[i].force.x);
00831         printf("%f, ", fmsg->slowForce[i].force.y);
00832         printf("%f>\n", fmsg->slowForce[i].force.z);
00833       }
00834       */
00835     }
00836 
00837     //CkPrintf("GBISENERGY[%i] c = % e, b = % e\n",timestep,coulEnergy, gbEnergy);
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   //printf("GBIS distributing forces COMPLETE!\n");
00860 }
00861 
00862 void ComputeGBISserMgr::recvForce(GBISForceMsg *msg) {
00863   //printf("GBIS recvForce()\n");
00864   gbisCompute->saveResults(msg);
00865   delete oldmsg;
00866   oldmsg = msg;
00867 }
00868 
00869 void ComputeGBISser::saveResults(GBISForceMsg *msg) {
00870   //printf("GBIS saveResults()\n");
00871   ResizeArrayIter<PatchElem> ap(patchList);
00872 
00873   ExtForce *results_ptr = msg->force;
00874   ExtForce *results_ptr_slow = msg->slowForce;
00875 
00876   // add in forces
00877   for (ap = ap.begin(); ap != ap.end(); ap++) {
00878     //CkPrintf("GBIS%i: ComputeGBIS(%i)::saveResults() openedForceBox",CkMyPe(),cid);
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     //CkPrintf("GBIS%i: slow[%i] = % e\n",CkMyPe(),i,sf[i].x);
00888       ++results_ptr;
00889       ++results_ptr_slow;
00890     }
00891     //CkPrintf("GBIS%i: ComputeGBIS(%i)::saveResults() closedForceBox",CkMyPe(),cid);
00892     (*ap).forceBox->close(&r);
00893   }
00894 
00895     //reduction->item(REDUCTION_ELECT_ENERGY) += msg->coulEnergy;
00896     reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbInterEnergy;
00897     reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbSelfEnergy;
00898     //CkPrintf("energies= % e, % e, % e\n",msg->coulEnergy, msg->gbInterEnergy, msg->gbSelfEnergy);
00899     reduction->submit();
00900 }
00901 
00902 /******************************************************************************
00903  * Send data to node 0
00904  ******************************************************************************/
00905 void ComputeGBISser::doWork()
00906 {
00907   //printf("GBIS doWork()\n");
00908   ResizeArrayIter<PatchElem> ap(patchList);
00909 
00910 #if 1
00911  // Skip computations if nothing to do.
00912  if ( ! patchList[0].p->flags.doNonbonded )
00913  {
00914    for (ap = ap.begin(); ap != ap.end(); ap++) {
00915      CompAtom *x = (*ap).positionBox->open();
00916      //CkPrintf("GBIS%i: ComputeGBIS(%i)::doWork() openedForceBox",CkMyPe(),cid);
00917      Results *r = (*ap).forceBox->open();
00918      (*ap).positionBox->close(&x);
00919      //CkPrintf("GBIS%i: ComputeGBIS(%i)::doWork() closedForceBox",CkMyPe(),cid);
00920      (*ap).forceBox->close(&r);
00921    }
00922    reduction->submit();
00923    return;
00924  }
00925 #endif
00926 
00927 
00928   // allocate message
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   //CkPrintf("SERIAL SLOW %i\n",msg->doSlow);
00941   msg->sequence = sequence();
00942 
00943   // get positions
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   Piecewise screening functions Hij dHij/drij
00978   r   distance
00979   r2  square distance
00980   ri  inverse distance
00981   rc  cutoff
00982   r0  descreened atom radius
00983   rs  descreening atom radius
00984   h   return value
00985   dh  return value
00986 */
00987 inline void h0Ser ( BigReal r, BigReal r2, BigReal ri,// 5.3%
00988 Real rc, BigReal r0, BigReal rs, BigReal & h ) {
00989   h = 0;
00990 }
00991 inline void dh0Ser ( BigReal r, BigReal r2, BigReal ri,// 5.3%
00992 Real rc, BigReal r0, BigReal rs, BigReal & dh ) {
00993   dh = 0;
00994 }
00995 
00996 inline void h1Ser ( BigReal r, BigReal r2, BigReal ri, //18.4%
00997 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
00998 
00999   BigReal rci = 1.0/rc;
01000   BigReal rmrs = r-rs;// 4 times
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, //18.4%
01009 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01010 
01011   BigReal rci = 1.0/rc;
01012   BigReal rmrs = r-rs;// 4 times
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,// 74.5%
01022 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01023 
01024     BigReal k = rs*ri; k*=k;//k=(rs/r)^2
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,// 74.5%
01028 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
01029 
01030     BigReal k = rs*ri; k*=k;//k=(rs/r)^2
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,// 1.4%
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,// 1.4%
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,// 0.4%
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,// 0.4%
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,// 0%, r<0.7Ang
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,// 0%, r<0.7Ang
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,//0%, one atom within other
01094 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
01095   h = 0;
01096 }
01097 inline void dh6Ser ( BigReal r, BigReal r2, BigReal ri,//0%, one atom within other
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 ) {//II
01105     h2Ser(r,r2,ri,rc,r0,rs,h); d = 2;
01106   } else if (r <= rc + rs && r > rc - rs) {//I
01107     h1Ser(r,r2,ri,rc,r0,rs,h); d = 1;
01108   } else if (r > rc + rs) {//0
01109     h0Ser(r,r2,ri,rc,r0,rs,h); d = 0;
01110   } else if( r <= 4*rs && r > r0 + rs ) {//III
01111     h3Ser(r,r2,ri,rc,r0,rs,h); d = 3;
01112   } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {//IV
01113     h4Ser(r,r2,ri,rc,r0,rs,h); d = 4;
01114   } else if (r0 < rs ) {//V
01115     h5Ser(r,r2,ri,rc,r0,rs,h); d = 5;
01116   } else {//VI
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 ) {//II
01123     dh2Ser(r,r2,ri,rc,r0,rs,dh); d = 2;
01124   } else if (r <= rc + rs && r > rc - rs) {//I
01125     dh1Ser(r,r2,ri,rc,r0,rs,dh); d = 1;
01126   } else if (r > rc + rs) {//0
01127     dh0Ser(r,r2,ri,rc,r0,rs,dh); d = 0;
01128   } else if( r <= 4*rs && r > r0 + rs ) {//III
01129     dh3Ser(r,r2,ri,rc,r0,rs,dh); d = 3;
01130   } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {//IV
01131     dh4Ser(r,r2,ri,rc,r0,rs,dh); d = 4;
01132   } else if (r0 < rs ) {//V
01133     dh5Ser(r,r2,ri,rc,r0,rs,dh); d = 5;
01134   } else {//VI
01135     dh6Ser(r,r2,ri,rc,r0,rs,dh); d = 6;
01136   }
01137 }
01138 inline void CalcHPairSer (
01139   BigReal r,//distance
01140   BigReal r2,//distance squared
01141   BigReal ri,//inverse distance
01142   BigReal rc,//cutoff
01143   BigReal ri0,
01144   BigReal rjs,
01145   BigReal rj0,
01146   BigReal ris,
01147   int & dij,//domain 1
01148   int & dji,//domain 2
01149   BigReal & hij,//output
01150   BigReal & hji//output
01151 ) {
01152   CalcHSer(r,r2,ri,rc,ri0,rjs,hij,dij);//hij
01153   CalcHSer(r,r2,ri,rc,rj0,ris,hji,dji);//hji
01154 }
01155 inline void CalcDHPairSer (
01156   BigReal r,//distance
01157   BigReal r2,
01158   BigReal ri,
01159   BigReal rc,//cutoff
01160   BigReal ri0,
01161   BigReal rjs,
01162   BigReal rj0,
01163   BigReal ris,
01164   int & dij,//domain 1
01165   int & dji,//domain 2
01166   BigReal & dhij,
01167   BigReal & dhji
01168 ) {
01169   CalcDHSer(r,r2,ri,rc,ri0,rjs,dhij,dij);//hij
01170   CalcDHSer(r,r2,ri,rc,rj0,ris,dhji,dji);//hji
01171 }
01172 
01173 
01174 /*******************************************************************************
01175 ********************************************************************************
01176 ***********  Phase 2 Inner Loop   **********************************************
01177 ********************************************************************************
01178 *******************************************************************************/
01179 
01180 /*
01181  * Calculate coulomb energy and force for single pair of atoms
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) {//not excluded
01193     //calculate Coulomb Energy
01194     coulE = -qiqj*epsilon_p_i*r_i;
01195 
01196     //calculate Coulomb Force
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  * Calculate GB Energy, GB dEdr force
01208  * also output intermediate values used in dEda
01209  */
01210 inline void Calc_dEdr_PairSer(//no longer does i==j
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,   //return
01226   BigReal & ddrGbE //return
01227 ) {
01228   //allocate local variables
01229   BigReal aiaj4,ddrDij,ddrf_i,ddrfij;
01230 
01231   //calculate GB energy
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;//dielectric term
01242   gbE = qiqj*Dij*f_i;
01243 
01244   //calculate energy derivatives
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  * Calculate summation element of dEda array
01253  * must calculate dEdr previously to retreive intermediate values
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,//return
01269   BigReal & dEdaj //return
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;//0
01275   dEdai = tmp_dEda/ai;
01276   dEdaj = tmp_dEda/aj;
01277 }
01278 
01279 /*
01280  * Calculate Coulomb and GB interaction and dEda element
01281  * for a pair of atoms
01282  */
01283 inline void Phase2_PairSer(//doesn't do self energies
01284 
01285 //input values
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 //return values
01300   BigReal & coulEij,
01301   BigReal & ddrCoulEij,
01302   BigReal & gbEij,
01303   BigReal & ddrGbEij,
01304   BigReal & dEdai,
01305   BigReal & dEdaj
01306 ) {
01307 
01308   //calculate Coulomb energy and force
01309   //Calc_Coul_Pair(r_i,qiqj,epsilon_p_i,exclij,scale14,coulEij,ddrCoulEij);
01310   coulEij = 0;
01311   ddrCoulEij = 0;
01312 
01313   //calculate GB energy and force
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   //calculate dEda
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