NAMD
ComputeGBISser.C
Go to the documentation of this file.
1 
7 /*******************************************************************************
8  *******************************************************************************
9  This serial version of GBIS is out of date and hasn't been vetted;
10  it is to be used as a way of testing new implicit solvent models
11  since all atomic coordinates are gathered into a single Compute.
12  *******************************************************************************
13  ******************************************************************************/
14 
15 /*
16  print debug
17  1 - once per step
18  2 - once per atom
19  3 - once per pair
20 */
21 //#define PRINT_SERFORCES
22 #define BENCH_PERIOD 1000
23 
24 //sum all forces
25 #define GBIS_DEDR_FORCE 1
26 #define GBIS_DEDA_FORCE 1
27 #define GBIS_COUL_FORCE 1
28 
29 #include "Vector.h"
30 #include <limits>
31 #include "InfoStream.h"
32 #include "Node.h"
33 #include "PatchMap.h"
34 #include "PatchMap.inl"
35 #include "AtomMap.h"
36 #include "ComputeGBISser.h"
37 #include "ComputeGBISserMgr.decl.h"
38 #include "PatchMgr.h"
39 #include "Molecule.h"
40 #include "ReductionMgr.h"
41 #include "ComputeMgr.h"
42 #include "ComputeMgr.decl.h"
43 #include "ComputeGBIS.inl"
44 // #define DEBUGM
45 #define MIN_DEBUG_LEVEL 3
46 #include "Debug.h"
47 #include "SimParameters.h"
48 #include "WorkDistrib.h"
49 #include "varsizemsg.h"
50 #include <stdlib.h>
51 #include <stdio.h>
52 #include <errno.h>
53 #include <time.h>
54 
55  struct vect {
56  BigReal x, y, z;
57  };
58 #ifndef COUL_CONST
59 #define COUL_CONST 332.0636 // ke [kcal*Ang/e^2]
60 #endif
61 
62 #ifndef TA
63 #define TA 0.333333333333333 // 1/3
64 #define TB 0.4 // 2/5
65 #define TC 0.428571428571428 // 3/7
66 #define TD 0.444444444444444 // 4/9
67 #define TE 0.454545454545454 // 5/11
68 #define DA 1.333333333333333 // 4* 1/3
69 #define DB 2.4 // 6* 2/5
70 #define DC 3.428571428571428 // 8* 3/7
71 #define DD 4.444444444444444 // 10*4/9
72 #define DE 5.454545454545454 // 12*5/11
73 #endif
74 
75 inline void Phase2_PairSer(//doesn't do self energies
76 
77 //input values
78  BigReal r,
79  BigReal r2,
80  BigReal r_i,
81  BigReal qiqj,
82  BigReal ai,
83  BigReal aj,
84  BigReal epsilon_p_i,
85  BigReal epsilon_s_i,
86  BigReal kappa,
87  int exclij,
88  BigReal scale14,
89  int stat,
90 
91 //return values
92  BigReal & coulEij,
93  BigReal & ddrCoulEij,
94  BigReal & gbEij,
95  BigReal & ddrGbEij,
96  BigReal & dEdai,
97  BigReal & dEdaj
98 );
99 
100 inline void CalcHPairSer (
101  BigReal r,//distance
102  BigReal r2,
103  BigReal ri,
104  BigReal rc,//cutoff
105  BigReal ri0,
106  BigReal rjs,
107  BigReal rj0,
108  BigReal ris,
109  int & dij,//domain 1
110  int & dji,//domain 2
111  BigReal & dhij,
112  BigReal & dhji);
113 inline void CalcDHPairSer (
114  BigReal r,//distance
115  BigReal r2,
116  BigReal ri,
117  BigReal rc,//cutoff
118  BigReal ri0,
119  BigReal rjs,
120  BigReal rj0,
121  BigReal ris,
122  int & dij,//domain 1
123  int & dji,//domain 2
124  BigReal & dhij,
125  BigReal & dhji);
126 
127 
131  Real rho;//coulomb radius
132  Real rho0;//coulomb radius
133  Real rhos;//coulomb radius
135  //int type;//atom type for S scaling table
136  int id;
137  int vdwType;
138 };
139 
140 class GBISCoordMsg : public CMessage_GBISCoordMsg {
141 public:
143  int numAtoms;
144  int doSlow;
146  int sequence;
147 };
148 
149 class GBISForceMsg : public CMessage_GBISForceMsg {
150 public:
156 };
157 
158 class ComputeGBISserMgr : public CBase_ComputeGBISserMgr {
159 public:
162 
163  void setCompute(ComputeGBISser *c) { gbisCompute = c; }
164  void recvCoord(GBISCoordMsg *);
165  void recvForce(GBISForceMsg *);
166 
167 private:
168  CProxy_ComputeGBISserMgr gbisProxy;
169  ComputeGBISser *gbisCompute;
170 
171  int numSources;
172  int numArrived;
173  GBISCoordMsg **coordMsgs;
174  int numAtoms;
175  ComputeGBISAtom *coord;
176  ExtForce *force;
177  ExtForce *slowForce;
178  GBISForceMsg *oldmsg;
179 
180  BigReal gbSelfEnergy;
181  BigReal gbInterEnergy;
182  BigReal coulEnergy;
183  int timestep;
184  clock_t t_start;
185  clock_t t_stop;
186  clock_t t1, t2;
187  double totalnamdtime;
188  double totalgbistime;
189  double inittime;
190  double psitime;
191  double alphatime;
192  double dEdrtime;
193  double dEdasumtime;
194  double dEdaprefixtime;
195  double dEdalooptime;
196  void calcGBISReg(int stat);
197  int blockSize;
198  int numBlocks;
199  int loop1Iter, loop2Iter, loop3Iter;
200  bool all2all;
201  BigReal MinErr;
202  BigReal MaxErr;
203  int numBins;
204  BigReal base;
205 };
206 
208  gbisProxy(thisgroup), gbisCompute(0), numSources(0), numArrived(0),
209  inittime(0),
210  psitime(0),
211  alphatime(0),
212  dEdrtime(0),
213  dEdasumtime(0),
214  dEdaprefixtime(0),
215  dEdalooptime(0), slowForce(0),
216  coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0), timestep(0) {
217  CkpvAccess(BOCclass_group).computeGBISserMgr = thisgroup;
218  t_start = clock();
219  t_stop = clock();
220  all2all = false;
221 }
222 
224  for ( int i=0; i<numSources; ++i ) { delete coordMsgs[i]; }
225  delete [] coordMsgs;
226  delete [] coord;
227  delete [] force;
228  delete [] slowForce;
229  delete oldmsg;
230 }
231 
234 {
235  CProxy_ComputeGBISserMgr::ckLocalBranch(
236  CkpvAccess(BOCclass_group).computeGBISserMgr)->setCompute(this);
237 
239 
240 }
241 
243 {
244 }
245 
246 /******************************************************************************
247  * scale from 1 to 0 between transition and cutoff
248  * s = 1 for [0,trans]
249  * s = s(r) for [trans, cutoff]
250  * s = 0 for [cutoff, infinity]
251  * meets the following 4 conditions:
252  * (1) s(trans) = 1
253  * (3) s'(trans) = 0
254  * (2) s(cutoff) = 0
255  * (4) s'(cutoff) = 0
256  ******************************************************************************/
257 inline void CalcScaleSer (
258  BigReal r,//rij
259  BigReal t,//transition distance
260  BigReal c,//cutoff distance
261  BigReal & s,//output s(r)
262  BigReal & d ) {//output s'(r)
263 
264  if (r <= t) { //[0,trans]
265  s = 1;
266  d = 0;
267  } else if (r < c) { //[trans,cutoff]
268  //precompute
269  BigReal ct = (c-t);
270  BigReal ct_i = 1/ct;
271  BigReal ct_i2 = ct_i*ct_i;
272  //BigReal ct_i4 = ct_i2*ct_i2;
273 
274  BigReal rt = r - t;
275  BigReal rt2 = rt*rt;
276  BigReal omrtct2 = 1-rt2*ct_i2;
277  s = omrtct2*omrtct2;
278 
279  BigReal rc = r - c;
280  d=s*4*rt/(rc*(rt+ct));
281 
282  } else { //[cutoff,infinity]
283  s = 0;
284  d = 0;
285  }
286 }
287 
288 
289 /**********************************************************
290  *
291  * Regular 3 sets of nested loops
292  *
293  * supports 1-2-4 scheme
294  *
295  **********************************************************/
296 void ComputeGBISserMgr::calcGBISReg(int stat) {
297  t1 = clock();
298  //printf("GBIS: ComputeGBIS::serial\n");
299 
300  /*********************************************************
301  * Initializations
302  ********************************************************/
303  Molecule *molecule = Node::Object()->molecule;
305  BigReal epsilon_s = simParams->solvent_dielectric;
306  BigReal epsilon_p = simParams->dielectric;
307  BigReal epsilon_s_i = 1/epsilon_s;
308  BigReal epsilon_p_i = 1/epsilon_p;
309  BigReal rho_0 = simParams->coulomb_radius_offset;
310  BigReal cutoff = simParams->cutoff;// + 1.8*0.96;//+max screen radius
311  BigReal trans = simParams->switchingDist;
312  BigReal fsMax = 1.728;
313  BigReal a_cut = simParams->alpha_cutoff-fsMax;//+max screen radius; partial sphere
314  //a_cut = 10;
315  BigReal a_cut2 = a_cut*a_cut;
316  BigReal a_cut_ps = a_cut + fsMax;//+max screen radius; partial sphere
317  BigReal cutoff2 = cutoff*cutoff;
318  BigReal a_cut_ps2 = a_cut_ps*a_cut_ps;
319  BigReal delta = simParams->gbis_delta;
320  BigReal beta = simParams->gbis_beta;
321  BigReal gamma = simParams->gbis_gamma;
322  BigReal alphaMax = simParams->alpha_max;
323  int exclude = simParams->exclude;
324  BigReal scale14 = (exclude==SCALED14) ? simParams->scale14 : 1.0;
325 
326 
327 #if PRINT > 0
328  //printf("cut=%e, rgbmax=%e, fsmax=%e, gbalpha=%e, gbbeta=%e, gbgamma=%e\n", cutoff, alphaMax, fsMax, delta, beta, gamma);
329 #endif
330 
331  //calc kappa
332  BigReal kappa = simParams->kappa;
333 
334  //debugging variables
335  gbInterEnergy = 0;
336  gbSelfEnergy = 0;
337  coulEnergy = 0;
338  BigReal fx, fy, fz;
339  loop1Iter = 0;
340  loop2Iter = 0;
341  loop3Iter = 0;
342  int loop1flops = 0;
343  int loop2flops = 0;
344  int loop3flops = 0;
345  BigReal loop1time = 0;
346  BigReal loop2time = 0;
347  BigReal loop3time = 0;
348 
349  Position ri, rj;
350  BigReal rhoi, rhoi0, rhois;
351  BigReal rhoj, rhoj0, rhojs;
352  BigReal dx, dy, dz;
353  BigReal r, r2, r_i, qiqj;
354  BigReal rnx, rny, rnz;
355  BigReal bornRad;
356  BigReal aiaj, aiaj4, expr2aiaj4, fij, finv, expkappa, Dij;
357  BigReal rc, shift, ddrshift;
358  BigReal tmp_dEda, dEdai, dEdaj, gbEij;
359  BigReal ddrfij, ddrfinv, ddrDij, ddrGbEij;
360  BigReal coulEij, ddrCoulEij;
361  BigReal forceCoul, forcedEdr, forceAlpha;
362  BigReal hij, hji, ddrHij, ddrHji;
363  BigReal tanhi, daidr, dajdr;
364  BigReal nbetapsi, gammapsi2;
365  vect *coulForce = new vect[numAtoms];
366  vect *gbForceR = new vect[numAtoms];//dEdr
367  vect *gbForceA = new vect[numAtoms];//dEda*dadr
368  BigReal *a = new BigReal[numAtoms];
369  BigReal *psi = new BigReal[numAtoms];
370  BigReal *dEda = new BigReal[numAtoms];
371  BigReal *ddrHijPrefix = new BigReal[numAtoms];
372 
373  //init arrays
374  for (int i = 0; i < numAtoms; i++) {
375  force[i].force.x =0.0;
376  force[i].force.y =0.0;
377  force[i].force.z =0.0;
378  slowForce[i].force.x =0.0;
379  slowForce[i].force.y =0.0;
380  slowForce[i].force.z =0.0;
381  dEda[i] = 0;
382 
383  psi[i] = 0;
384  ddrHijPrefix[i] = 0;
385  coulForce[i].x = 0;
386  coulForce[i].y = 0;
387  coulForce[i].z = 0;
388  gbForceR[i].x = 0;//dEdr
389  gbForceR[i].y = 0;//dEdr
390  gbForceR[i].z = 0;//dEdr
391  gbForceA[i].x = 0;//dEda*dadr
392  gbForceA[i].y = 0;//dEda*dadr
393  gbForceA[i].z = 0;//dEda*dadr
394  }
395 
396  t2 = clock();
397  inittime = (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
398 
399  /**********************************************************
400  *
401  * Loop 1 : accumulate Hij into psi array
402  * cutoff = a_cut = 25A
403  * every 2 timesteps
404  *
405  **********************************************************/
406  double missTime = 0;
407  double hitTime = 0;
408  t1 = clock();
409  int dij, dji;//which hij domain 0123456 - for debugging only
410  BigReal dhij, dhji;
411  //int TESTITER = 1;
412  //for (int k = 0; k < TESTITER; k++)
413  for (int i = 0; i < numAtoms; i++) {
414  ri = coord[i].position;
415  for (int j = i+1; j < numAtoms; j++) {
416  rj = coord[j].position;
417  dx = (ri.x - rj.x);
418  dy = (ri.y - rj.y);
419  dz = (ri.z - rj.z);
420  r2 = dx*dx+dy*dy+dz*dz;
421  if (r2 > a_cut_ps2) continue;
422  r_i = 1.0/sqrt(r2);;
423  r = 1.0/r_i;;
424  loop1Iter++;
425  BigReal rhoi0 = coord[i].rho0;
426  BigReal rhojs = coord[j].rhos;
427  BigReal rhoj0 = coord[j].rho0;
428  BigReal rhois = coord[i].rhos;
429  CalcHPairSer(r,r2,r_i,a_cut, coord[i].rho0, coord[j].rhos,
430  coord[j].rho0, coord[i].rhos,dij,dji,hij,hji);
431 
432 #ifdef PRINT_COMP
433  CkPrintf("PSI(%04i)[%04i,%04i] = %i%i % .4e % .4e\n",timestep,coord[i].id,coord[j].id,dij,dji,hij,hji);
434  //CkPrintf("S_Hij_%i\n",dij);
435  //CkPrintf("S_Hji_%i\n",dji);
436 #endif
437  psi[i] += hij;
438  psi[j] += hji;
439 
440  }
441  }
442  t2 = clock();
443  psitime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
444  loop1time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
445  t1 = clock();
446 
447 
448 
449  /**********************************************************
450  *
451  * Atom 1: Calculate alpha based on phi for each atom
452  *
453  **********************************************************/
454  BigReal totPsi = 0;
455  for (int i = 0; i < numAtoms; i++) {
456  //CkPrintf("Srho[%i] = %.3e %.3e\n",coord[i].id,coord[i].rho0,coord[i].rhos);
457 
458  rhoi0 = coord[i].rho0;
459  rhoi = coord[i].rho;
460  //CkPrintf("GBIS_SER: psi[%03i] = %.4e\n",i,psi[i]);
461  totPsi += psi[i];
462  psi[i] *= rhoi0;
463  bornRad=1/(1/rhoi0-1/rhoi*tanh(psi[i]*(delta+psi[i]*(-beta+gamma*psi[i]))));
464  bornRad = (1.0/bornRad < 1.0/alphaMax) ? alphaMax : bornRad;
465  a[i] = bornRad;
466 #ifdef PRINT_COMP
467  CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",timestep,coord[i].id,bornRad);
468 #endif
469  }
470  t2 = clock();
471  alphatime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
472  t1 = clock();
473 
474 
475 
476  /**********************************************************
477  *
478  * Loop 2 : dEda
479  *
480  **********************************************************/
481  for (int i = 0; i < numAtoms; i++) {
482  ri = coord[i].position;
483  for (int j = i+1; j < numAtoms; j++) {
484  rj = coord[j].position;
485  dx = (ri.x - rj.x);//rptI
486  dy = (ri.y - rj.y);//rptI
487  dz = (ri.z - rj.z);//rptI
488  r2 = dx*dx+dy*dy+dz*dz;//rptI
489  if (r2 > cutoff2) continue;
490 
491  //calculate distance
492  loop2Iter++;
493  qiqj = (-COUL_CONST*coord[i].charge);
494  qiqj *= coord[j].charge;
495  r_i = 1.0/sqrt(r2);
496  r = 1.0/r_i;//rptI
497  rnx = dx*r_i;//normalized vector
498  rny = dy*r_i;
499  rnz = dz*r_i;
500 
501  //calculate scaling for energy smoothing
502  BigReal scale, ddrScale;
503  //CalcScaleSer(r, trans, cutoff, scale, ddrScale);
504  if ( false ) {
505  BigReal ratio = r2 / cutoff2;
506  scale = ratio - 1;
507  scale *= scale;
508  ddrScale = 4.0*r*(r2-cutoff2)/cutoff2/cutoff2;
509  } else {
510  scale = 1;
511  ddrScale = 0;
512  }
513 
514 
515  //BigReal ai = a[i];
516  //BigReal aj = a[j];
517  int exclij = 0;//molecule->checkexcl(i,j);
518  //CkPrintf("GBIS_SER_excl[%i,%i] %i\n",i,j,exclij);
519  //CkPrintf("GBIS_DOFULL = %i\n",stat);
520 
521  Phase2_PairSer(r,r2,r_i,qiqj,a[i],a[j],epsilon_p_i,epsilon_s_i,kappa,exclij,
522  scale14,stat,coulEij,ddrCoulEij,gbEij,ddrGbEij,dEdai,dEdaj);
523  int id1 = coord[i].id;
524  int id2 = coord[j].id;
525  if (id1 > id2 ) {
526  int tmp = id2;
527  id2 = id1;
528  id1 = tmp;
529  }//print ids as lower index, higher index
530 
531  //accumulate energies
532  gbInterEnergy += gbEij *scale;
533  forcedEdr = -(ddrGbEij)*scale-(gbEij)*ddrScale;
534 
535  //add GB force
536  fx = rnx*forcedEdr;
537  fy = rny*forcedEdr;
538  fz = rnz*forcedEdr;
539 #ifdef PRINT_COMP
540  CkPrintf("DEDR(%04i)[%04i,%04i] = % .4e\n",timestep,i,j,forcedEdr);
541  CkPrintf("DASM(%04i)[%04i,%04i] = % .4e % .4e\n",timestep,i,j,dEdai*scale,dEdaj*scale);
542  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);
543 
544 #endif
545  gbForceR[i].x += fx;
546  gbForceR[i].y += fy;
547  gbForceR[i].z += fz;
548  gbForceR[j].x -= fx;
549  gbForceR[j].y -= fy;
550  gbForceR[j].z -= fz;
551 
552  //add dEda
553  if (stat) {
554  dEda[i] += dEdai*scale;
555  dEda[j] += dEdaj*scale;
556  }
557 
558  }//end j loop
559  }//end i loop
560  t2 = clock();
561  dEdasumtime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
562  loop2time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
563  t1 = clock();
564 
565 
566 
567 /*******************************************************************************
568 *
569 * Atom 2 : Calculate dHij Prefix
570 *
571 *******************************************************************************/
572  for (int i = 0; i < numAtoms; i++) {
573  //add diagonal dEda term
574  fij = a[i];//inf
575  expkappa = exp(-kappa*fij);//0
576  Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
577  //add diagonal GB energy term
578  gbEij = -COUL_CONST*coord[i].charge*coord[i].charge*Dij/fij;
579  gbSelfEnergy += 0.5*gbEij;//self energy
580 //CkPrintf("self_energy[%03i] = %e\n",i,0.5*gbEij);
581  //CkPrintf("gbSelfEnergy[%03i] = % e\n", coord[i].id,0.5*gbEij);
582 
583  //calculate dHij prefix
584  if (stat) {
585  dEdai = -0.5*COUL_CONST*coord[i].charge*coord[i].charge
586  *(kappa*epsilon_s_i*expkappa-Dij/fij)/a[i];
587  //CkPrintf("SER_dEdai[%03i]%i = % e\n",i,timestep,dEdai);
588  dEda[i] += dEdai;
589  BigReal dedasum = dEda[i];
590  //CkPrintf("SER_dEdaSum[%03i]%i = % e\n",i,timestep,dEda[i]);
591 
592  rhoi0 = coord[i].rho0;
593  rhoi = rhoi0+rho_0;
594  BigReal psii = psi[i];
595  nbetapsi = -beta*psii;
596  gammapsi2 = gamma*psii*psii;
597  tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
598  daidr = a[i]*a[i]*rhoi0/rhoi*(1-tanhi*tanhi)
599  * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
600  ddrHijPrefix[i] = daidr*dEda[i];
601 #ifdef PRINT_COMP
602  CkPrintf("DHDR(%04i)[%04i] = % .4e\n",timestep,coord[i].id, ddrHijPrefix[i]);
603 #endif
604  }
605  }
606  t2 = clock();
607  dEdaprefixtime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
608  //gbEnergy = gbInterEnergy + gbSelfEnergy;
609  t1 = clock();
610 
611  /**********************************************************
612  *
613  * Loop 3 : Calculate dEda*dadr Forces
614  *
615  **********************************************************/
616  if (stat) {
617  BigReal dhij, dhji;
618  int dij, dji;
619  //for (int k = 0; k < TESTITER; k++)
620  for (int i = 0; i < numAtoms; i++) {
621  ri = coord[i].position;
622  //CkPrintf("SER_dHdrPrefix[%i]%i = % .5e\n",coord[i].id,timestep,ddrHijPrefix[i]);
623 
624  for (int j = i+1; j < numAtoms; j++) {
625  rj = coord[j].position;
626  dx = (ri.x - rj.x);//rptI
627  dy = (ri.y - rj.y);//rptI
628  dz = (ri.z - rj.z);//rptI
629  r2 = dx*dx+dy*dy+dz*dz;//rptI
630  //loop3flops += 9;
631  if (r2 > a_cut_ps2) continue;
632  r_i = 1.0/sqrt(r2);//rptI
633  loop3Iter++;
634  //calculate dHij for pair
635  r = 1.0/r_i;
636  CalcDHPairSer(r,r2,r_i,a_cut,
637  coord[i].rho0,coord[j].rhos,
638  coord[j].rho0,coord[i].rhos,
639  dij,dji,dhij,dhji);
640 
641  //calculate and add dEijdai,j*dai,jdrij force
642  forceAlpha = -r_i*(ddrHijPrefix[i]*dhij+ddrHijPrefix[j]*dhji);
643 #ifdef PRINT_COMP
644  CkPrintf("DEDA(%04i)[%04i,%04i] = %i%i % .4e % .4e % .4e\n",timestep,i,j,dij,dji,dhij,dhji, forceAlpha/r_i);
645 #endif
646  fx = dx * forceAlpha;
647  fy = dy * forceAlpha;
648  fz = dz * forceAlpha;
649 
650  gbForceA[i].x += fx;
651  gbForceA[i].y += fy;
652  gbForceA[i].z += fz;
653  gbForceA[j].x -= fx;
654  gbForceA[j].y -= fy;
655  gbForceA[j].z -= fz;
656  }//end inner summation loop
657  }//end outer summation loop
658  }//end if stat
659  t2 = clock();
660  dEdalooptime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
661  loop3time += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
662  t1 = clock();
663 
664  //add forces
665  for (int i = 0; i < numAtoms; i++) {
666  force[i].force.x =0.0;
667  force[i].force.y =0.0;
668  force[i].force.z =0.0;
669  slowForce[i].force.x =0.0;
670  slowForce[i].force.y =0.0;
671  slowForce[i].force.z =0.0;
672 
673 #if GBIS_COUL_FORCE
674  force[i].force.x += coulForce[i].x;
675  force[i].force.y += coulForce[i].y;
676  force[i].force.z += coulForce[i].z;
677 #endif
678 
679 #if GBIS_DEDR_FORCE
680  force[i].force.x += gbForceR[i].x;
681  force[i].force.y += gbForceR[i].y;
682  force[i].force.z += gbForceR[i].z;
683 #endif
684 
685 #if GBIS_DEDA_FORCE
686  if (stat) {
687  slowForce[i].force.x += gbForceA[i].x;
688  slowForce[i].force.y += gbForceA[i].y;
689  slowForce[i].force.z += gbForceA[i].z;
690  //CkPrintf("SERIAL SLOW %e %e %e\n",gbForceA[i].x,gbForceA[i].y,gbForceA[i].z);
691  }
692  //force[i].force.x += gbForceA[i].x;
693  //force[i].force.y += gbForceA[i].y;
694  //force[i].force.z += gbForceA[i].z;
695 #endif
696 
697 
698 #ifdef PRINT_SERFORCES
699  BigReal fC, fR, fA;
700  fx = 0;
701  fy = 0;
702  fz = 0;
703 #if GBIS_COUL_FORCE
704  fx += coulForce[i].x;
705  fy += coulForce[i].y;
706  fz += coulForce[i].z;
707  //fC = sqrt(fx*fx+fy*fy+fz*fz);
708  //fprintf(stderr, "%i %e %e %e ",i,fx, fy, fz);
709 #endif
710 #if GBIS_DEDR_FORCE
711  fx += gbForceR[i].x;
712  fy += gbForceR[i].y;
713  fz += gbForceR[i].z;
714  //fR = sqrt(fx*fx+fy*fy+fz*fz);
715  //fprintf(stderr, "%e %e %e",fx, fy, fz);
716 #endif
717 #if GBIS_DEDA_FORCE
718  fx += gbForceA[i].x;
719  fy += gbForceA[i].y;
720  fz += gbForceA[i].z;
721  //fA = sqrt(fx*fx+fy*fy+fz*fz);
722  //fprintf(stderr, "%5i % .5e % .5e % .5e\n",i, fx, fy, fz);
723 #endif
724  //fprintf(stderr, "%5i % .5e % .5e % .5e\n",i, fx, fy, fz);
725 #endif //if print>1
726  }
727  t2 = clock();
728  inittime += (double)((double)t2-(double)t1)/CLOCKS_PER_SEC;
729 
730  timestep++;
731 
732  delete[] coulForce;
733  delete[] gbForceR;
734  delete[] gbForceA;
735  delete[] a;
736  delete[] psi;
737  delete[] dEda;
738  delete[] ddrHijPrefix;
739  //printf("GBIS calcReg COMPLETE!\n");
740 }
741 
743  //printf("GBIS recvCoord()\n");
744  if ( ! numSources ) {//init receiving coord
745  numSources = (PatchMap::Object())->numNodesWithPatches();
746  coordMsgs = new GBISCoordMsg*[numSources];
747  for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
748  numArrived = 0;
749  numAtoms = Node::Object()->molecule->numAtoms;
750  coord = new ComputeGBISAtom[numAtoms];
751  force = new ExtForce[numAtoms];
752  slowForce = new ExtForce[numAtoms];
753  timestep = msg[0].sequence;
754  }
755 
756 
757  //receive coord
758  int i;
759  for ( i=0; i < msg->numAtoms; ++i ) {
760  coord[msg->coord[i].id] = msg->coord[i];
761  }
762 
763  coordMsgs[numArrived] = msg;
764  ++numArrived;
765 
766  if ( numArrived < numSources ) return;
767  numArrived = 0;
768 
769 
770  /**********************************************************
771  *
772  * All sources arrived; calculate energy, forces
773  *
774  **********************************************************/
775  t_start = clock();
776  //how long did namd take to cycle back here
777  double namdtime = (double)((double)t_start-(double)t_stop)/CLOCKS_PER_SEC;
778  totalnamdtime += namdtime;
779  //choice choose
780  calcGBISReg(msg->doSlow);
781 
782  t_stop = clock();
783  double gbistime = (double)((double)t_stop-(double)t_start)/CLOCKS_PER_SEC;
784  //printf("GBIS: elapsednamd(%i)=%f\n", timestep, namdtime);
785  //printf("GBIS: elapsedgbis(%i)=%f\n", timestep, gbistime);
786  totalgbistime += gbistime;
787  //printf("GBIS: total(%i)=%f\n", timestep, totaltime);
788  if (timestep % BENCH_PERIOD == 0) {
789  printf("\n");
790  printf("GBIS: t_GB=%f sec for %i steps\n",totalgbistime,BENCH_PERIOD);
791  printf("GBIS: t_MD=%f sec for %i steps\n",totalnamdtime,BENCH_PERIOD);
792  printf("GBIS: init=%f sec for %i steps\n", inittime, BENCH_PERIOD);
793  printf("GBIS: psi=%f sec for %i steps\n", psitime, BENCH_PERIOD);
794  printf("GBIS: alpha=%f sec for %i steps\n", alphatime, BENCH_PERIOD);
795  printf("GBIS: dEdasum=%f sec for %i steps\n", dEdasumtime, BENCH_PERIOD);
796  printf("GBIS: dEdaprefix=%f sec for %i steps\n",dEdaprefixtime,BENCH_PERIOD);
797  printf("GBIS: dEdaloop=%f sec for %i steps\n", dEdalooptime,BENCH_PERIOD);
798  printf("GBIS: loop1=%i iters\n", loop1Iter);
799  printf("GBIS: loop2=%i iters\n", loop2Iter);
800  printf("GBIS: loop3=%i iters\n", loop3Iter);
801  printf("\n");
802  totalgbistime = 0;
803  totalnamdtime = 0;
804  inittime = 0;
805  psitime = 0;
806  alphatime = 0;
807  dEdrtime = 0;
808  dEdasumtime = 0;
809  dEdaprefixtime = 0;
810  dEdalooptime = 0;
811  }
812 
813  // distribute forces
814  //printf("GBIS distributing forces\n");
815  for ( int j=0; j < numSources; ++j ) {
816  GBISCoordMsg *cmsg = coordMsgs[j];
817  coordMsgs[j] = 0;
818  GBISForceMsg *fmsg = new (cmsg->numAtoms, cmsg->numAtoms, 0) GBISForceMsg;
819  //does this init slowForce, force?
820  for ( int i=0; i < cmsg->numAtoms; ++i ) {
821  fmsg->force[i] = force[cmsg->coord[i].id];
822  fmsg->slowForce[i] = slowForce[cmsg->coord[i].id];
823  /*
824  if (i == 1000 ) {
825  printf("%i force: <", i);
826  printf("%f, ", fmsg->force[i].force.x);
827  printf("%f, ", fmsg->force[i].force.y);
828  printf("%f>\n", fmsg->force[i].force.z);
829  printf("%i slwfr: <", i);
830  printf("%f, ", fmsg->slowForce[i].force.x);
831  printf("%f, ", fmsg->slowForce[i].force.y);
832  printf("%f>\n", fmsg->slowForce[i].force.z);
833  }
834  */
835  }
836 
837  //CkPrintf("GBISENERGY[%i] c = % e, b = % e\n",timestep,coulEnergy, gbEnergy);
838  if ( ! j ) {
839 #if GBIS_DEDR_FORCE
840  fmsg->gbSelfEnergy = gbSelfEnergy;
841  fmsg->gbInterEnergy = gbInterEnergy;
842 #else
843  fmsg->gbSelfEnergy = 0;
844  fmsg->gbInterEnergy = 0;
845 #endif
846 #if GBIS_COUL_FORCE
847  fmsg->coulEnergy = coulEnergy;
848 #else
849  fmsg->coulEnergy = 0;
850 #endif
851  } else {
852  fmsg->gbSelfEnergy = 0;
853  fmsg->gbInterEnergy = 0;
854  fmsg->coulEnergy = 0;
855  }
856  gbisProxy[cmsg->sourceNode].recvForce(fmsg);
857  delete cmsg;
858  }
859  //printf("GBIS distributing forces COMPLETE!\n");
860 }
861 
863  //printf("GBIS recvForce()\n");
864  gbisCompute->saveResults(msg);
865  delete oldmsg;
866  oldmsg = msg;
867 }
868 
870  //printf("GBIS saveResults()\n");
872 
873  ExtForce *results_ptr = msg->force;
874  ExtForce *results_ptr_slow = msg->slowForce;
875 
876  // add in forces
877  for (ap = ap.begin(); ap != ap.end(); ap++) {
878  //CkPrintf("GBIS%i: ComputeGBIS(%i)::saveResults() openedForceBox",CkMyPe(),cid);
879  Results *r = (*ap).forceBox->open();
880  Force *f = r->f[Results::nbond];
881  Force *sf = r->f[Results::slow];
882  int numAtoms = (*ap).p->getNumAtoms();
883 
884  for(int i=0; i<numAtoms; ++i) {
885  f[i] += results_ptr->force;
886  sf[i] += results_ptr_slow->force;
887  //CkPrintf("GBIS%i: slow[%i] = % e\n",CkMyPe(),i,sf[i].x);
888  ++results_ptr;
889  ++results_ptr_slow;
890  }
891  //CkPrintf("GBIS%i: ComputeGBIS(%i)::saveResults() closedForceBox",CkMyPe(),cid);
892  (*ap).forceBox->close(&r);
893  }
894 
895  //reduction->item(REDUCTION_ELECT_ENERGY) += msg->coulEnergy;
896  reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbInterEnergy;
897  reduction->item(REDUCTION_ELECT_ENERGY) += msg->gbSelfEnergy;
898  //CkPrintf("energies= % e, % e, % e\n",msg->coulEnergy, msg->gbInterEnergy, msg->gbSelfEnergy);
899  reduction->submit();
900 }
901 
902 /******************************************************************************
903  * Send data to node 0
904  ******************************************************************************/
906 {
907  //printf("GBIS doWork()\n");
909 
910 #if 1
911  // Skip computations if nothing to do.
912  if ( ! patchList[0].p->flags.doNonbonded )
913  {
914  for (ap = ap.begin(); ap != ap.end(); ap++) {
915  CompAtom *x = (*ap).positionBox->open();
916  //CkPrintf("GBIS%i: ComputeGBIS(%i)::doWork() openedForceBox",CkMyPe(),cid);
917  Results *r = (*ap).forceBox->open();
918  (*ap).positionBox->close(&x);
919  //CkPrintf("GBIS%i: ComputeGBIS(%i)::doWork() closedForceBox",CkMyPe(),cid);
920  (*ap).forceBox->close(&r);
921  }
922  reduction->submit();
923  return;
924  }
925 #endif
926 
927 
928  // allocate message
929  int numLocalAtoms = 0;
930  for (ap = ap.begin(); ap != ap.end(); ap++) {
931  numLocalAtoms += (*ap).p->getNumAtoms();
932  }
933 
934  GBISCoordMsg *msg = new (numLocalAtoms, 0) GBISCoordMsg;
935  msg->sourceNode = CkMyPe();
936  msg->numAtoms = numLocalAtoms;
937  ComputeGBISAtom *data_ptr = msg->coord;
938  SimParameters *simParams = Node::Object()->simParameters;
939  msg->doSlow = patchList[0].p->flags.doFullElectrostatics;
940  //CkPrintf("SERIAL SLOW %i\n",msg->doSlow);
941  msg->sequence = sequence();
942 
943  // get positions
944  for (ap = ap.begin(); ap != ap.end(); ap++) {
945  FullAtomList &atoms = (*ap).p->getAtomList();
946  CompAtom *x = (*ap).positionBox->open();
947  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
948  int numAtoms = (*ap).p->getNumAtoms();
949  for(int i=0; i<numAtoms; ++i)
950  {
951  data_ptr->position = x[i].position;
952  data_ptr->charge = x[i].charge;
953  data_ptr->mass = atoms[i].mass;
954  data_ptr->id = xExt[i].id;
955  data_ptr->rho = MassToRadius(data_ptr->mass);
956  SimParameters *simParams = Node::Object()->simParameters;
957  data_ptr->rho0 = data_ptr->rho - simParams->coulomb_radius_offset;
958  data_ptr->rhos = data_ptr->rho0 * MassToScreen(data_ptr->mass);
959  data_ptr->vdwType = x[i].vdwType;
960  ++data_ptr;
961  }
962 
963 #if 0
964  if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
965  else { (*ap).positionBox->close(&x); }
966 #endif
967  (*ap).positionBox->close(&x);
968  }
969 
970  CProxy_ComputeGBISserMgr gbisProxy(CkpvAccess(BOCclass_group).computeGBISserMgr);
971  gbisProxy[0].recvCoord(msg);
972 
973 }
974 #include "ComputeGBISserMgr.def.h"
975 
976 /*
977  Piecewise screening functions Hij dHij/drij
978  r distance
979  r2 square distance
980  ri inverse distance
981  rc cutoff
982  r0 descreened atom radius
983  rs descreening atom radius
984  h return value
985  dh return value
986 */
987 inline void h0Ser ( BigReal r, BigReal r2, BigReal ri,// 5.3%
988 Real rc, BigReal r0, BigReal rs, BigReal & h ) {
989  h = 0;
990 }
991 inline void dh0Ser ( BigReal r, BigReal r2, BigReal ri,// 5.3%
992 Real rc, BigReal r0, BigReal rs, BigReal & dh ) {
993  dh = 0;
994 }
995 
996 inline void h1Ser ( BigReal r, BigReal r2, BigReal ri, //18.4%
997 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
998 
999  BigReal rci = 1.0/rc;
1000  BigReal rmrs = r-rs;// 4 times
1001  BigReal rmrsi = 1.0/rmrs;
1002  BigReal rmrs2 = rmrs*rmrs;
1003  BigReal rs2 = rs*rs;
1004  BigReal logr = log(rmrs*rci);
1005  BigReal rci2 = rci*rci;
1006  h = 0.125*ri*(1 + 2*r*rmrsi + rci2*(r2 - 4*rc*r - rs2) + 2*logr);
1007 }
1008 inline void dh1Ser ( BigReal r, BigReal r2, BigReal ri, //18.4%
1009 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1010 
1011  BigReal rci = 1.0/rc;
1012  BigReal rmrs = r-rs;// 4 times
1013  BigReal rmrsi = 1.0/rmrs;
1014  BigReal rmrs2 = rmrs*rmrs;
1015  BigReal rs2 = rs*rs;
1016  BigReal logr = log(rmrs*rci);
1017  BigReal rci2 = rci*rci;
1018  dh = ri*ri*(-0.25*logr - (rc*rc - rmrs2)*(rs2 + r2)*0.125*rci2*rmrsi*rmrsi);
1019 }
1020 
1021 inline void h2Ser ( BigReal r, BigReal r2, BigReal ri,// 74.5%
1022 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
1023 
1024  BigReal k = rs*ri; k*=k;//k=(rs/r)^2
1025  h = rs*ri*ri*k*(TA+k*(TB+k*(TC+k*(TD+k*TE))));
1026 }
1027 inline void dh2Ser ( BigReal r, BigReal r2, BigReal ri,// 74.5%
1028 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1029 
1030  BigReal k = rs*ri; k*=k;//k=(rs/r)^2
1031  dh = -rs*ri*ri*ri*k*(DA+k*(DB+k*(DC+k*(DD+k*DE))));
1032 }
1033 
1034 inline void h3Ser ( BigReal r, BigReal r2, BigReal ri,// 1.4%
1035 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
1036  BigReal r2mrs2i = 1.0/(r2-rs*rs);
1037  h = 0.5 * ( rs*r2mrs2i + 0.5 * log((r-rs)/(r+rs))*ri );
1038 }
1039 inline void dh3Ser ( BigReal r, BigReal r2, BigReal ri,// 1.4%
1040 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1041  BigReal rs2 = rs*rs;
1042  BigReal r2mrs2i = 1.0/(r2-rs2);
1043  dh = -0.25*ri*(2*(r2+rs2)*rs*r2mrs2i*r2mrs2i + ri*log((r-rs)/(r+rs)));
1044 }
1045 
1046 inline void h4Ser ( BigReal r, BigReal r2, BigReal ri,// 0.4%
1047 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
1048  BigReal ri2 = ri*ri;
1049  BigReal r02 = r0*r0;
1050  BigReal rs2 = rs*rs;
1051  BigReal r0i = 1.0/r0;
1052  BigReal rspri = 1.0/(r+rs);
1053  BigReal logr = log(r0*rspri);
1054  BigReal r02mrs2 = r02-rs2;
1055  BigReal rilogr = ri*logr;
1056  h = 0.25*( r0i*(2-0.5*(r0i*ri*(r2 + r02 - rs2))) - rspri + rilogr );
1057 }
1058 inline void dh4Ser ( BigReal r, BigReal r2, BigReal ri,// 0.4%
1059 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1060  BigReal ri2 = ri*ri;
1061  BigReal r02 = r0*r0;
1062  BigReal rs2 = rs*rs;
1063  BigReal r0i = 1.0/r0;
1064  BigReal rspri = 1.0/(r+rs);
1065  BigReal logr = log(r0*rspri);
1066  BigReal r02mrs2 = r02-rs2;
1067  BigReal rilogr = ri*logr;
1068  dh = 0.25*( (-0.5+(r2*r02mrs2 - 2*r*rs*rs2+rs2*r02mrs2)
1069  * 0.5*ri2*rspri*rspri)*r0i*r0i - ri*rilogr );
1070 }
1071 
1072 inline void h5Ser ( BigReal r, BigReal r2, BigReal ri,// 0%, r<0.7Ang
1073 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
1074  BigReal rs2 = rs*rs;
1075  BigReal r2mrs2i = 1/(r2-rs2);
1076  BigReal rsr2mrs2i = rs*r2mrs2i;
1077  BigReal rprs = r+rs;
1078  BigReal rmrs = r-rs;
1079  BigReal logr = 0.5*ri*log(-rmrs/rprs);
1080  h = 0.5*( rsr2mrs2i + 2/r0 + logr );
1081 }
1082 inline void dh5Ser ( BigReal r, BigReal r2, BigReal ri,// 0%, r<0.7Ang
1083 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1084  BigReal rs2 = rs*rs;
1085  BigReal r2mrs2i = 1/(r2-rs2);
1086  BigReal rsr2mrs2i = rs*r2mrs2i;
1087  BigReal rprs = r+rs;
1088  BigReal rmrs = r-rs;
1089  BigReal logr = 0.5*ri*log(-rmrs/rprs);
1090  dh = -0.5*ri*((rs2+r2)*rsr2mrs2i*r2mrs2i+logr );
1091 }
1092 
1093 inline void h6Ser ( BigReal r, BigReal r2, BigReal ri,//0%, one atom within other
1094 BigReal rc, BigReal r0, BigReal rs, BigReal & h ) {
1095  h = 0;
1096 }
1097 inline void dh6Ser ( BigReal r, BigReal r2, BigReal ri,//0%, one atom within other
1098 BigReal rc, BigReal r0, BigReal rs, BigReal & dh ) {
1099  dh = 0;
1100 }
1101 
1102 inline void CalcHSer ( BigReal r, BigReal r2, BigReal ri,
1103 BigReal rc, BigReal r0, BigReal rs, BigReal & h, int & d) {
1104  if( r <= rc - rs && r > 4*rs ) {//II
1105  h2Ser(r,r2,ri,rc,r0,rs,h); d = 2;
1106  } else if (r <= rc + rs && r > rc - rs) {//I
1107  h1Ser(r,r2,ri,rc,r0,rs,h); d = 1;
1108  } else if (r > rc + rs) {//0
1109  h0Ser(r,r2,ri,rc,r0,rs,h); d = 0;
1110  } else if( r <= 4*rs && r > r0 + rs ) {//III
1111  h3Ser(r,r2,ri,rc,r0,rs,h); d = 3;
1112  } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {//IV
1113  h4Ser(r,r2,ri,rc,r0,rs,h); d = 4;
1114  } else if (r0 < rs ) {//V
1115  h5Ser(r,r2,ri,rc,r0,rs,h); d = 5;
1116  } else {//VI
1117  h6Ser(r,r2,ri,rc,r0,rs,h); d = 6;
1118  }
1119 }
1120 inline void CalcDHSer ( BigReal r, BigReal r2, BigReal ri,
1121 BigReal rc, BigReal r0, BigReal rs, BigReal & dh, int & d) {
1122  if( r <= rc - rs && r > 4*rs ) {//II
1123  dh2Ser(r,r2,ri,rc,r0,rs,dh); d = 2;
1124  } else if (r <= rc + rs && r > rc - rs) {//I
1125  dh1Ser(r,r2,ri,rc,r0,rs,dh); d = 1;
1126  } else if (r > rc + rs) {//0
1127  dh0Ser(r,r2,ri,rc,r0,rs,dh); d = 0;
1128  } else if( r <= 4*rs && r > r0 + rs ) {//III
1129  dh3Ser(r,r2,ri,rc,r0,rs,dh); d = 3;
1130  } else if ( r <= r0 + rs && r > (r0>rs?r0-rs:rs-r0) ) {//IV
1131  dh4Ser(r,r2,ri,rc,r0,rs,dh); d = 4;
1132  } else if (r0 < rs ) {//V
1133  dh5Ser(r,r2,ri,rc,r0,rs,dh); d = 5;
1134  } else {//VI
1135  dh6Ser(r,r2,ri,rc,r0,rs,dh); d = 6;
1136  }
1137 }
1138 inline void CalcHPairSer (
1139  BigReal r,//distance
1140  BigReal r2,//distance squared
1141  BigReal ri,//inverse distance
1142  BigReal rc,//cutoff
1143  BigReal ri0,
1144  BigReal rjs,
1145  BigReal rj0,
1146  BigReal ris,
1147  int & dij,//domain 1
1148  int & dji,//domain 2
1149  BigReal & hij,//output
1150  BigReal & hji//output
1151 ) {
1152  CalcHSer(r,r2,ri,rc,ri0,rjs,hij,dij);//hij
1153  CalcHSer(r,r2,ri,rc,rj0,ris,hji,dji);//hji
1154 }
1155 inline void CalcDHPairSer (
1156  BigReal r,//distance
1157  BigReal r2,
1158  BigReal ri,
1159  BigReal rc,//cutoff
1160  BigReal ri0,
1161  BigReal rjs,
1162  BigReal rj0,
1163  BigReal ris,
1164  int & dij,//domain 1
1165  int & dji,//domain 2
1166  BigReal & dhij,
1167  BigReal & dhji
1168 ) {
1169  CalcDHSer(r,r2,ri,rc,ri0,rjs,dhij,dij);//hij
1170  CalcDHSer(r,r2,ri,rc,rj0,ris,dhji,dji);//hji
1171 }
1172 
1173 
1174 /*******************************************************************************
1175 ********************************************************************************
1176 *********** Phase 2 Inner Loop **********************************************
1177 ********************************************************************************
1178 *******************************************************************************/
1179 
1180 /*
1181  * Calculate coulomb energy and force for single pair of atoms
1182  */
1183 inline void Calc_Coul_PairSer(
1184  BigReal r_i,
1185  BigReal qiqj,
1186  BigReal epsilon_p_i,
1187  int exclij,
1188  BigReal scale14,
1189  BigReal & coulE,
1190  BigReal & ddrCoulE
1191 ) {
1192  if (exclij != EXCHCK_FULL) {//not excluded
1193  //calculate Coulomb Energy
1194  coulE = -qiqj*epsilon_p_i*r_i;
1195 
1196  //calculate Coulomb Force
1197  if (exclij == EXCHCK_MOD)
1198  coulE *= scale14;
1199  ddrCoulE = -r_i*coulE;
1200  } else {
1201  coulE = 0;
1202  ddrCoulE = 0;
1203  }
1204 }
1205 
1206 /*
1207  * Calculate GB Energy, GB dEdr force
1208  * also output intermediate values used in dEda
1209  */
1210 inline void Calc_dEdr_PairSer(//no longer does i==j
1211  BigReal r,
1212  BigReal r2,
1213  BigReal qiqj,
1214  BigReal ai,
1215  BigReal aj,
1216  BigReal kappa,
1217  BigReal epsilon_p_i,
1218  BigReal epsilon_s_i,
1219  BigReal & aiaj,
1220  BigReal & expr2aiaj4,
1221  BigReal & fij,
1222  BigReal & f_i,
1223  BigReal & expkappa,
1224  BigReal & Dij,
1225  BigReal & gbE, //return
1226  BigReal & ddrGbE //return
1227 ) {
1228  //allocate local variables
1229  BigReal aiaj4,ddrDij,ddrf_i,ddrfij;
1230 
1231  //calculate GB energy
1232  aiaj = ai*aj;
1233  aiaj4 = 4*aiaj;
1234  expr2aiaj4 = exp(-r2/aiaj4);
1235  fij = sqrt(r2+aiaj*expr2aiaj4);
1236  f_i = 1/fij;
1237  if (kappa > 0)
1238  expkappa = exp(-kappa*fij);
1239  else
1240  expkappa = 1.0;
1241  Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
1242  gbE = qiqj*Dij*f_i;
1243 
1244  //calculate energy derivatives
1245  ddrfij = r*f_i*(1 - 0.25*expr2aiaj4);
1246  ddrf_i = -ddrfij*f_i*f_i;
1247  ddrDij = kappa*expkappa*ddrfij*epsilon_s_i;
1248  ddrGbE = qiqj*(ddrDij*f_i+Dij*ddrf_i);
1249 }
1250 
1251 /*
1252  * Calculate summation element of dEda array
1253  * must calculate dEdr previously to retreive intermediate values
1254  */
1255 inline void Calc_dEda_PairSer(
1256  BigReal r2,
1257  BigReal ai,
1258  BigReal aj,
1259  BigReal qiqj,
1260  BigReal kappa,
1261  BigReal aiaj,
1262  BigReal expkappa,
1263  BigReal expr2aiaj4,
1264  BigReal fij,
1265  BigReal f_i,
1266  BigReal Dij,
1267  BigReal epsilon_s_i,
1268  BigReal & dEdai,//return
1269  BigReal & dEdaj //return
1270 ) {
1271 
1272  BigReal tmp_dEda = 0.5*qiqj*f_i*f_i
1273  *(kappa*epsilon_s_i*expkappa-Dij*f_i)
1274  *(aiaj+0.25*r2)*expr2aiaj4;//0
1275  dEdai = tmp_dEda/ai;
1276  dEdaj = tmp_dEda/aj;
1277 }
1278 
1279 /*
1280  * Calculate Coulomb and GB interaction and dEda element
1281  * for a pair of atoms
1282  */
1283 inline void Phase2_PairSer(//doesn't do self energies
1284 
1285 //input values
1286  BigReal r,
1287  BigReal r2,
1288  BigReal r_i,
1289  BigReal qiqj,
1290  BigReal ai,
1291  BigReal aj,
1292  BigReal epsilon_p_i,
1293  BigReal epsilon_s_i,
1294  BigReal kappa,
1295  int exclij,
1296  BigReal scale14,
1297  int doSlow,
1298 
1299 //return values
1300  BigReal & coulEij,
1301  BigReal & ddrCoulEij,
1302  BigReal & gbEij,
1303  BigReal & ddrGbEij,
1304  BigReal & dEdai,
1305  BigReal & dEdaj
1306 ) {
1307 
1308  //calculate Coulomb energy and force
1309  //Calc_Coul_Pair(r_i,qiqj,epsilon_p_i,exclij,scale14,coulEij,ddrCoulEij);
1310  coulEij = 0;
1311  ddrCoulEij = 0;
1312 
1313  //calculate GB energy and force
1314  BigReal aiaj,expr2aiaj4,fij,f_i,expkappa,Dij;
1315  Calc_dEdr_PairSer(r,r2,qiqj,ai,aj,kappa,epsilon_p_i,epsilon_s_i,
1316  aiaj,expr2aiaj4,fij,f_i,expkappa,Dij,gbEij,ddrGbEij);
1317 
1318  //calculate dEda
1319  if (doSlow) {
1320  Calc_dEda_PairSer(r2,ai,aj,qiqj,kappa,aiaj,expkappa,expr2aiaj4,
1321  fij,f_i,Dij,epsilon_s_i,dEdai,dEdaj);
1322  }
1323 }
1324 
static Node * Object()
Definition: Node.h:86
#define TE
Definition: ComputeGBIS.inl:36
#define SCALED14
Definition: SimParameters.h:43
void h2Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
virtual ~ComputeGBISser()
int sequence(void)
Definition: Compute.h:64
#define TD
Definition: ComputeGBIS.inl:35
void saveResults(GBISForceMsg *)
BigReal solvent_dielectric
#define TA
Definition: ComputeGBIS.inl:32
int ComputeID
Definition: NamdTypes.h:183
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
static __thread atom * atoms
float Real
Definition: common.h:107
BigReal & item(int i)
Definition: ReductionMgr.h:312
void Phase2_PairSer(BigReal r, BigReal r2, BigReal r_i, BigReal qiqj, BigReal ai, BigReal aj, BigReal epsilon_p_i, BigReal epsilon_s_i, BigReal kappa, int exclij, BigReal scale14, int stat, BigReal &coulEij, BigReal &ddrCoulEij, BigReal &gbEij, BigReal &ddrGbEij, BigReal &dEdai, BigReal &dEdaj)
BigReal gbis_gamma
BigReal z
Definition: Vector.h:66
void dh3Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
ComputeGBISser(ComputeID c)
Position position
Definition: NamdTypes.h:53
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
void h3Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
#define TC
Definition: ComputeGBIS.inl:34
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
void recvForce(GBISForceMsg *)
void Calc_dEdr_PairSer(BigReal r, BigReal r2, BigReal qiqj, BigReal ai, BigReal aj, BigReal kappa, BigReal epsilon_p_i, BigReal epsilon_s_i, BigReal &aiaj, BigReal &expr2aiaj4, BigReal &fij, BigReal &f_i, BigReal &expkappa, BigReal &Dij, BigReal &gbE, BigReal &ddrGbE)
#define DA
Definition: ComputeGBIS.inl:37
#define EXCHCK_MOD
Definition: Molecule.h:86
void dh6Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
BigReal switchingDist
BigReal coulomb_radius_offset
BigReal y
#define DD
Definition: ComputeGBIS.inl:40
#define COUL_CONST
BigReal x
void dh2Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
Charge charge
Definition: NamdTypes.h:54
void dh4Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
void h5Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
void h0Ser(BigReal r, BigReal r2, BigReal ri, Real rc, BigReal r0, BigReal rs, BigReal &h)
BigReal alpha_max
ResizeArrayIter< T > end(void) const
void recvCoord(GBISCoordMsg *)
void dh0Ser(BigReal r, BigReal r2, BigReal ri, Real rc, BigReal r0, BigReal rs, BigReal &dh)
void dh1Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
void Calc_Coul_PairSer(BigReal r_i, BigReal qiqj, BigReal epsilon_p_i, int exclij, BigReal scale14, BigReal &coulE, BigReal &ddrCoulE)
BigReal gbSelfEnergy
void Calc_dEda_PairSer(BigReal r2, BigReal ai, BigReal aj, BigReal qiqj, BigReal kappa, BigReal aiaj, BigReal expkappa, BigReal expr2aiaj4, BigReal fij, BigReal f_i, BigReal Dij, BigReal epsilon_s_i, BigReal &dEdai, BigReal &dEdaj)
void dh5Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh)
#define DB
Definition: ComputeGBIS.inl:38
#define EXCHCK_FULL
Definition: Molecule.h:85
void CalcHSer(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h, int &d)
void CalcDHSer(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &dh, int &d)
#define DC
Definition: ComputeGBIS.inl:39
Force * f[maxNumForces]
Definition: PatchTypes.h:67
BigReal coulEnergy
BigReal gbis_beta
void CalcHPairSer(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal ri0, BigReal rjs, BigReal rj0, BigReal ris, int &dij, int &dji, BigReal &dhij, BigReal &dhji)
BigReal x
Definition: Vector.h:66
BigReal alpha_cutoff
int numAtoms
Definition: Molecule.h:556
Force force
Definition: NamdTypes.h:202
#define simParams
Definition: Output.C:127
short vdwType
Definition: NamdTypes.h:55
ExtForce * slowForce
BigReal gbInterEnergy
BigReal y
Definition: Vector.h:66
static float MassToRadius(Mass mi)
Definition: ComputeGBIS.inl:55
static float MassToScreen(Mass mi)
void h4Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
float Mass
Definition: ComputeGBIS.inl:20
void h1Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
BigReal gbis_delta
BigReal dielectric
void submit(void)
Definition: ReductionMgr.h:323
ComputeGBISAtom * coord
#define TB
Definition: ComputeGBIS.inl:33
void h6Ser(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal r0, BigReal rs, BigReal &h)
void setCompute(ComputeGBISser *c)
gridSize x
ExclusionSettings exclude
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cutoff2
Molecule * molecule
Definition: Node.h:176
#define BENCH_PERIOD
ExtForce * force
#define DE
Definition: ComputeGBIS.inl:41
ResizeArrayIter< T > begin(void) const
BigReal z
double BigReal
Definition: common.h:112
void CalcDHPairSer(BigReal r, BigReal r2, BigReal ri, BigReal rc, BigReal ri0, BigReal rjs, BigReal rj0, BigReal ris, int &dij, int &dji, BigReal &dhij, BigReal &dhji)
for(int i=0;i< n1;++i)
void CalcScaleSer(BigReal r, BigReal t, BigReal c, BigReal &s, BigReal &d)