NAMD
ComputeNonbondedUtil.C
Go to the documentation of this file.
1 
7 /*
8  Common operations for ComputeNonbonded classes
9 */
10 
11 #ifdef WIN32
12 extern "C" {
13  double erfc(double);
14 }
15 #endif
16 
17 #include "InfoStream.h"
18 #include "ComputeNonbondedUtil.h"
19 #include "SimParameters.h"
20 #include "Node.h"
21 #include "Molecule.h"
22 #include "LJTable.h"
23 #include "ReductionMgr.h"
24 #include "Parameters.h"
25 #include "MsmMacros.h"
26 #include <stdio.h>
27 
28 #ifdef NAMD_CUDA
30 #endif
31 
32 #ifdef NAMD_MIC
33  extern void send_build_mic_force_table();
34 #endif
35 
63 #if defined(NAMD_MIC)
64  BigReal* ComputeNonbondedUtil::mic_table_base_ptr;
65  int ComputeNonbondedUtil::mic_table_n;
66  int ComputeNonbondedUtil::mic_table_n_16;
67 #endif
68 #ifdef NAMD_AVXTILES
69 int ComputeNonbondedUtil::avxTilesMode;
70 float* ComputeNonbondedUtil::avx_tiles_eps4_sigma = 0;
71 float* ComputeNonbondedUtil::avx_tiles_eps4_sigma_14 = 0;
72 #endif
73 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
74 float* ComputeNonbondedUtil::knl_table_alloc;
75 float* ComputeNonbondedUtil::knl_fast_ener_table;
76 float* ComputeNonbondedUtil::knl_fast_grad_table;
77 float* ComputeNonbondedUtil::knl_scor_ener_table;
78 float* ComputeNonbondedUtil::knl_scor_grad_table;
79 float* ComputeNonbondedUtil::knl_slow_ener_table;
80 float* ComputeNonbondedUtil::knl_slow_grad_table;
81 float* ComputeNonbondedUtil::knl_excl_ener_table;
82 float* ComputeNonbondedUtil::knl_excl_grad_table;
83 #ifdef NAMD_KNL
84 float* ComputeNonbondedUtil::knl_corr_ener_table;
85 float* ComputeNonbondedUtil::knl_corr_grad_table;
86 #endif
87 #endif
119 // BigReal ComputeNonbondedUtil::d0;
120 // fepb
127 //fepe
131 
133 
136 
142 
144 
146 
149 
151 
152 // Ported by JLai -- JE - Go
155 int ComputeNonbondedUtil::goMethod; //6.3.11
156 // End of port -- JLai
157 
162 
167 
172 
177 
178 // define splitting function
179 #define SPLIT_NONE 1
180 #define SPLIT_SHIFT 2
181 #define SPLIT_C1 3
182 #define SPLIT_XPLOR 4
183 #define SPLIT_C2 5
184 #define SPLIT_MARTINI 6
185 
187 {
190  reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
192  reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
193  // Ported by JLai
194  reduction->item(REDUCTION_GRO_LJ_ENERGY) += data[groLJEnergyIndex];
198  // End of port -- JLai
199 //fepb
200  reduction->item(REDUCTION_ELECT_ENERGY_F) += data[electEnergyIndex_s];
202  reduction->item(REDUCTION_LJ_ENERGY_F) += data[vdwEnergyIndex_s];
203 
210 //fepe
211  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
212  ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,fullElectVirialIndex);
213  ADD_VECTOR(reduction,REDUCTION_PAIR_VDW_FORCE,data,pairVDWForceIndex);
214  ADD_VECTOR(reduction,REDUCTION_PAIR_ELECT_FORCE,data,pairElectForceIndex);
215  reduction->item(REDUCTION_COMPUTE_CHECKSUM) += 1.;
216 }
217 
219  SubmitReduction *reduction)
220 {
221  if (!reduction) return;
222  int numAtomTypes = pressureProfileAtomTypes;
223  // For ease of calculation we stored interactions between types
224  // i and j in (ni+j). For efficiency now we coalesce the
225  // cross interactions so that just i<=j are stored.
226  const int arraysize = 3*pressureProfileSlabs;
227  size_t nelems = arraysize*(numAtomTypes*(numAtomTypes+1))/2;
228  BigReal *arr = new BigReal[nelems];
229  memset(arr, 0, nelems*sizeof(BigReal));
230 
231  int i, j;
232  for (i=0; i<numAtomTypes; i++) {
233  for (j=0; j<numAtomTypes; j++) {
234  int ii=i;
235  int jj=j;
236  if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
237  const int reductionOffset = (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
238  for (int k=0; k<arraysize; k++) {
239  arr[reductionOffset+k] += data[k];
240  }
241  data += arraysize;
242  }
243  }
244  // copy into reduction
245  reduction->add(nelems, arr);
246  delete [] arr;
247 }
248 
250  NAMD_bug("Tried to call missing nonbonded compute routine.");
251 }
252 
254 {
255  if ( CkMyRank() ) return;
256 
257  // These defaults die cleanly if nothing appropriate is assigned.
274 
276  Parameters * params = Node::Object()->parameters;
277 
278  table_ener = params->table_ener;
279  rowsize = params->rowsize;
280  columnsize = params->columnsize;
281 
282  commOnly = simParams->commOnly;
283  fixedAtomsOn = ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces );
284 
285  qmForcesOn = simParams->qmForcesOn ;
286 
287  cutoff = simParams->cutoff;
289  cutoff2_f = cutoff2;
290 
291 //fepb
292  alchFepOn = simParams->alchFepOn;
293  alchThermIntOn = simParams->alchThermIntOn;
294  alchWCAOn = simParams->alchWCAOn;
295  alchDecouple = simParams->alchDecouple;
296 
297  lesOn = simParams->lesOn;
298  lesScaling = lesFactor = 0;
299 
300  Bool tabulatedEnergies = simParams->tabulatedEnergies;
303 
304  delete [] lambda_table;
305  lambda_table = 0;
306 
310 
311  // Ported by JLai -- Original JE - Go
312  goGroPair = simParams->goGroPair;
313  goForcesOn = simParams->goForcesOn;
314  goMethod = simParams->goMethod;
315  // End of port
316 
317  accelMDOn = simParams->accelMDOn;
318 
319  drudeNbthole = simParams->drudeOn && (simParams->drudeNbtholeCut > 0.0);
320 
321  if ( drudeNbthole ) {
322 #ifdef NAMD_CUDA
323  NAMD_die("drudeNbthole is not supported in CUDA version");
324 #endif
325  if ( lesOn )
326  NAMD_die("drudeNbthole is not supported with locally enhanced sampling");
327  if ( pairInteractionOn )
328  NAMD_die("drudeNbthole is not supported with pair interaction calculation");
329  if ( pressureProfileOn )
330  NAMD_die("drudeNbthole is not supported with pressure profile calculation");
331  }
332 
333  if ( alchFepOn ) {
334 #ifdef NAMD_CUDA
335  NAMD_die("Alchemical free-energy perturbation is not supported in CUDA version");
336 #endif
353  } else if ( alchThermIntOn ) {
354 #ifdef NAMD_CUDA
355  NAMD_die("Alchemical thermodynamic integration is not supported in CUDA version");
356 #endif
373  } else if ( lesOn ) {
374 #ifdef NAMD_CUDA
375  NAMD_die("Locally enhanced sampling is not supported in CUDA version");
376 #endif
377  lesFactor = simParams->lesFactor;
378  lesScaling = 1.0 / (double)lesFactor;
379  lambda_table = new BigReal[(lesFactor+1)*(lesFactor+1)];
380  for ( int ip=0; ip<=lesFactor; ++ip ) {
381  for ( int jp=0; jp<=lesFactor; ++jp ) {
382  BigReal lambda_pair = 1.0;
383  if (ip || jp ) {
384  if (ip && jp && ip != jp) {
385  lambda_pair = 0.0;
386  } else {
387  lambda_pair = lesScaling;
388  }
389  }
390  lambda_table[(lesFactor+1)*ip+jp] = lambda_pair;
391  }
392  }
409  } else if ( pressureProfileOn) {
410 #ifdef NAMD_CUDA
411  NAMD_die("Pressure profile calculation is not supported in CUDA version");
412 #endif
415 
432  } else if ( pairInteractionOn ) {
433 #ifdef NAMD_CUDA
434  NAMD_die("Pair interaction calculation is not supported in CUDA version");
435 #endif
442  } else if ( tabulatedEnergies ) {
443 #ifdef NAMD_CUDA
444  NAMD_die("Tabulated energies is not supported in CUDA version");
445 #endif
462  } else if ( goForcesOn ) {
463 #ifdef NAMD_CUDA
464  NAMD_die("Go forces is not supported in CUDA version");
465 #endif
482  } else {
499  }
500 
501 //fepe
502 
503  dielectric_1 = 1.0/simParams->dielectric;
504  if ( simParams->soluteScalingOn ) {
505  delete ljTable;
506  ljTable = 0;
507  }
508  if ( ! ljTable ) ljTable = new LJTable;
510  scaling = simParams->nonbondedScaling;
511  if ( simParams->exclude == SCALED14 )
512  {
513  scale14 = simParams->scale14;
514  }
515  else
516  {
517  scale14 = 1.;
518  }
519  if ( simParams->switchingActive )
520  {
521  switchOn = simParams->switchingDist;
522  switchOn_1 = 1.0/switchOn;
523  // d0 = 1.0/(cutoff-switchOn);
525  c0 = 1.0/(cutoff2-switchOn2);
526 
527  if ( simParams->vdwForceSwitching ) {
528  double switchOn3 = switchOn * switchOn2;
529  double cutoff3 = cutoff * cutoff2;
530  double switchOn6 = switchOn3 * switchOn3;
531  double cutoff6 = cutoff3 * cutoff3;
532  v_vdwa_f = v_vdwa = -1. / ( switchOn6 * cutoff6 );
533  v_vdwb_f = v_vdwb = -1. / ( switchOn3 * cutoff3 );
534  k_vdwa_f = k_vdwa = cutoff6 / ( cutoff6 - switchOn6 );
535  k_vdwb_f = k_vdwb = cutoff3 / ( cutoff3 - switchOn3 );
536  cutoff_3_f = cutoff_3 = 1. / cutoff3;
537  cutoff_6_f = cutoff_6 = 1. / cutoff6;
538 
539  } else if ( simParams->martiniSwitching ) { // switching fxn for Martini RBCG
540 
541  BigReal p6 = 6;
542  BigReal A6 = p6 * ((p6+1)*switchOn-(p6+4)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,2));
543  BigReal B6 = -p6 * ((p6+1)*switchOn-(p6+3)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,3));
544  BigReal C6 = 1.0/pow(cutoff,p6)-A6/3.0*pow(cutoff-switchOn,3)-B6/4.0*pow(cutoff-switchOn,4);
545 
546  BigReal p12 = 12;
547  BigReal A12 = p12 * ((p12+1)*switchOn-(p12+4)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,2));
548  BigReal B12 = -p12 * ((p12+1)*switchOn-(p12+3)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,3));
549  BigReal C12 = 1.0/pow(cutoff,p12)-A12/3.0*pow(cutoff-switchOn,3)-B12/4.0*pow(cutoff-switchOn,4);
550 
551  A6_f = A6; B6_f = B6; C6_f = C6;
552  A12_f = A12; B12_f = B12; C12_f = C12;
554 
555  }
556 
557  }
558  else
559  {
560  switchOn = cutoff;
561  switchOn_1 = 1.0/switchOn;
562  // d0 = 0.; // avoid division by zero
564  c0 = 0.; // avoid division by zero
565  }
566  c1 = c0*c0*c0;
567  c3 = 3.0 * (cutoff2 - switchOn2);
568  c5 = 0;
569  c6 = 0;
570  c7 = 0;
571  c8 = 0;
572 
573  const int PMEOn = simParams->PMEOn;
574  const int MSMOn = simParams->MSMOn;
575  const int MSMSplit = simParams->MSMSplit;
576 
577  if ( PMEOn ) {
578  ewaldcof = simParams->PMEEwaldCoefficient;
579  BigReal TwoBySqrtPi = 1.12837916709551;
580  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
581  }
582 
583  int splitType = SPLIT_NONE;
584  if ( simParams->switchingActive ) splitType = SPLIT_SHIFT;
585  if ( simParams->martiniSwitching ) splitType = SPLIT_MARTINI;
586  if ( simParams->fullDirectOn || simParams->FMAOn || PMEOn || MSMOn ||
587  simParams->FMMOn ) {
588  switch ( simParams->longSplitting ) {
589  case C2:
590  splitType = SPLIT_C2;
591  break;
592 
593  case C1:
594  splitType = SPLIT_C1;
595  break;
596 
597  case XPLOR:
598  NAMD_die("Sorry, XPLOR splitting not supported.");
599  break;
600 
601  case SHARP:
602  NAMD_die("Sorry, SHARP splitting not supported.");
603  break;
604 
605  default:
606  NAMD_die("Unknown splitting type found!");
607 
608  }
609  }
610 
611  BigReal r2_tol = 0.1;
612 
613  r2_delta = 1.0;
614  r2_delta_exp = 0;
615  while ( r2_delta > r2_tol ) { r2_delta /= 2.0; r2_delta_exp += 1; }
616  r2_delta_1 = 1.0 / r2_delta;
617 
618  if ( ! CkMyPe() ) {
619  iout << iINFO << "NONBONDED TABLE R-SQUARED SPACING: " <<
620  r2_delta << "\n" << endi;
621  }
622 
623  BigReal r2_tmp = 1.0;
624  int cutoff2_exp = 0;
625  while ( (cutoff2 + r2_delta) > r2_tmp ) { r2_tmp *= 2.0; cutoff2_exp += 1; }
626 
627  int i;
628  int n = (r2_delta_exp + cutoff2_exp) * 64 + 1;
629  table_length = n;
630  #if defined(NAMD_MIC)
631  int n_16 = (n + 15) & (~15);
632  #endif
633 
634  if ( ! CkMyPe() ) {
635  iout << iINFO << "NONBONDED TABLE SIZE: " <<
636  n << " POINTS\n" << endi;
637  }
638 
639  if ( table_alloc ) delete [] table_alloc;
640  #if defined(NAMD_MIC)
641  table_alloc = new BigReal[61*n_16+16];
642  BigReal *table_align = table_alloc;
643  while ( ((long)table_align) % 128 ) ++table_align;
644  mic_table_base_ptr = table_align;
645  mic_table_n = n;
646  mic_table_n_16 = n_16;
647  table_noshort = table_align;
648  table_short = table_align + 16*n_16;
649  slow_table = table_align + 32*n_16;
650  fast_table = table_align + 36*n_16;
651  scor_table = table_align + 40*n_16;
652  corr_table = table_align + 44*n_16;
653  full_table = table_align + 48*n_16;
654  vdwa_table = table_align + 52*n_16;
655  vdwb_table = table_align + 56*n_16;
656  r2_table = table_align + 60*n_16;
657  #else
658  table_alloc = new BigReal[61*n+16];
659  BigReal *table_align = table_alloc;
660  while ( ((long)table_align) % 128 ) ++table_align;
661  table_noshort = table_align;
662  table_short = table_align + 16*n;
663  slow_table = table_align + 32*n;
664  fast_table = table_align + 36*n;
665  scor_table = table_align + 40*n;
666  corr_table = table_align + 44*n;
667  full_table = table_align + 48*n;
668  vdwa_table = table_align + 52*n;
669  vdwb_table = table_align + 56*n;
670  r2_table = table_align + 60*n;
671  #endif
672  BigReal *fast_i = fast_table + 4;
673  BigReal *scor_i = scor_table + 4;
674  BigReal *slow_i = slow_table + 4;
675  BigReal *vdwa_i = vdwa_table + 4;
676  BigReal *vdwb_i = vdwb_table + 4;
677  BigReal *r2_i = r2_table; *(r2_i++) = r2_delta;
678  BigReal r2_limit = simParams->limitDist * simParams->limitDist;
679  if ( r2_limit < r2_delta ) r2_limit = r2_delta;
680  int r2_delta_i = 0; // entry for r2 == r2_delta
681 
682 #ifdef NAMD_AVXTILES
683  {
684  avxTilesMode = 1;
685  if (!simParams->useAVXTiles)
686  avxTilesMode = 0;
687  if (avxTilesMode) {
688  if (simParams->vdwGeometricSigma) avxTilesMode = 2;
689  if ( simParams->fullDirectOn || simParams->FMAOn || PMEOn || MSMOn ||
690  simParams->FMMOn ) {
691  if ( simParams->longSplitting != C1)
692  avxTilesMode = 3;
693  } else
694  avxTilesMode = 3;
695  if (avxTilesMode == 1) {
696  const int table_dim = ljTable->get_table_dim();
697  Real A, B, A14, B14;
698  for (int i = 0; i < table_dim; i++)
699  for (int j = i+1; j < table_dim; j++)
700  if (params->get_vdw_pair_params(i, j, &A, &B, &A14, &B14))
701  avxTilesMode = 2;
702  }
703  if (avxTilesMode > 1)
704  iout << iINFO << "AVX-512 TILES WILL USE SHORT-RANGE INTERPOLATION ("
705  << avxTilesMode << ")\n";
706  else {
707  Parameters *params = Node::Object()->parameters;
708  const int num_params = params->get_num_vdw_params();
709  if ( avx_tiles_eps4_sigma ) delete [] avx_tiles_eps4_sigma;
710  if ( avx_tiles_eps4_sigma_14 ) delete [] avx_tiles_eps4_sigma_14;
711  avx_tiles_eps4_sigma = new float[num_params*2];
712  avx_tiles_eps4_sigma_14 = new float[num_params*2];
713  for (int i = 0; i < num_params; i++) {
714  Real sigma, sigma_14, epsilon, epsilon_14;
715  params->get_vdw_params(&sigma, &epsilon, &sigma_14,
716  &epsilon_14,i);
717  avx_tiles_eps4_sigma[i*2] = 4.0 * scaling * epsilon;
718  avx_tiles_eps4_sigma[i*2 + 1] = sigma;
719  avx_tiles_eps4_sigma_14[i*2] = 4.0 * scaling * epsilon_14;
720  avx_tiles_eps4_sigma_14[i*2 + 1] = sigma_14;
721  }
722  }
723  }
724  }
725 #endif
726 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
727  if ( knl_table_alloc ) delete [] knl_table_alloc;
728  if ( KNL_TABLE_MAX_R_1 < 1.f || KNL_TABLE_FACTOR < 1 ||
729  KNL_TABLE_FACTOR !=
730  static_cast<int>(1.0 / KNL_TABLE_MAX_R_1 * KNL_TABLE_SIZE))
731  NAMD_bug("Inconsistent KNL preprocessor settings.");
732 #ifdef NAMD_KNL
733  knl_table_alloc = new float[10*KNL_TABLE_SIZE];
734 #else
735  knl_table_alloc = new float[8*KNL_TABLE_SIZE];
736 #endif
737  knl_fast_ener_table = knl_table_alloc;
738  knl_fast_grad_table = knl_table_alloc + KNL_TABLE_SIZE;
739  knl_scor_ener_table = knl_table_alloc + 2*KNL_TABLE_SIZE;
740  knl_scor_grad_table = knl_table_alloc + 3*KNL_TABLE_SIZE;
741  knl_slow_ener_table = knl_table_alloc + 4*KNL_TABLE_SIZE;
742  knl_slow_grad_table = knl_table_alloc + 5*KNL_TABLE_SIZE;
743  knl_excl_ener_table = knl_table_alloc + 6*KNL_TABLE_SIZE;
744  knl_excl_grad_table = knl_table_alloc + 7*KNL_TABLE_SIZE;
745  knl_fast_ener_table[0] = 0.;
746  knl_fast_grad_table[0] = 0.;
747  knl_scor_ener_table[0] = 0.;
748  knl_scor_grad_table[0] = 0.;
749  knl_slow_ener_table[0] = 0.;
750  knl_slow_grad_table[0] = 0.;
751  knl_excl_ener_table[0] = 0.;
752  knl_excl_grad_table[0] = 0.;
753 #ifdef NAMD_KNL
754  knl_corr_ener_table = knl_table_alloc + 8*KNL_TABLE_SIZE;
755  knl_corr_grad_table = knl_table_alloc + 9*KNL_TABLE_SIZE;
756  knl_corr_ener_table[0] = 0.;
757  knl_corr_grad_table[0] = 0.;
758 #endif
759  for ( int knl_table = 0; knl_table < 2; ++knl_table ) {
760  int nn = n;
761  if ( knl_table ) {
762  nn = KNL_TABLE_SIZE-1;
763  }
764  for ( i=1; i<nn; ++i ) {
765 #else
766  // fill in the table, fix up i==0 (r2==0) below
767  for ( i=1; i<n; ++i ) {
768 #endif
769 
770  const BigReal r2_base = r2_delta * ( 1 << (i/64) );
771  const BigReal r2_del = r2_base / 64.0;
772  BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
773 
774  BigReal r = sqrt(r2);
775 
776 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
777  if ( knl_table ) {
778  r = (double)(KNL_TABLE_FACTOR-2)/(double)(i);
779  r2 = r*r;
780  } else
781 #endif
782  if ( r2 <= r2_limit ) r2_delta_i = i;
783 
784  const BigReal r_1 = 1.0/r;
785  const BigReal r_2 = 1.0/r2;
786 
787  // fast_ is defined as (full_ - slow_)
788  // corr_ and fast_ are both zero at the cutoff, full_ is not
789  // all three are approx 1/r at short distances
790 
791  // for actual interpolation, we use fast_ for fast forces and
792  // scor_ = slow_ + corr_ - full_ and slow_ for slow forces
793  // since these last two are of small magnitude
794 
795  BigReal fast_energy, fast_gradient;
796  BigReal scor_energy, scor_gradient;
797  BigReal slow_energy, slow_gradient;
798 
799  // corr_ is PME direct sum, or similar correction term
800  // corr_energy is multiplied by r until later
801  // corr_gradient is multiplied by -r^2 until later
802  BigReal corr_energy, corr_gradient;
803 
804 
805  if ( PMEOn ) {
806  BigReal tmp_a = r * ewaldcof;
807  BigReal tmp_b = erfc(tmp_a);
808  corr_energy = tmp_b;
809  corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
810  } else if ( MSMOn ) {
811  BigReal a_1 = 1.0/cutoff;
812  BigReal r_a = r * a_1;
813  BigReal g, dg;
814  SPOLY(&g, &dg, r_a, MSMSplit);
815  corr_energy = 1 - r_a * g;
816  corr_gradient = 1 + r_a*r_a * dg;
817  } else {
818  corr_energy = corr_gradient = 0;
819  }
820 
821  switch(splitType) {
822  case SPLIT_NONE:
823  fast_energy = 1.0/r;
824  fast_gradient = -1.0/r2;
825  scor_energy = scor_gradient = 0;
826  slow_energy = slow_gradient = 0;
827  break;
828  case SPLIT_SHIFT: {
829  BigReal shiftVal = r2/cutoff2 - 1.0;
830  shiftVal *= shiftVal;
831  BigReal dShiftVal = 2.0 * (r2/cutoff2 - 1.0) * 2.0*r/cutoff2;
832  fast_energy = shiftVal/r;
833  fast_gradient = dShiftVal/r - shiftVal/r2;
834  scor_energy = scor_gradient = 0;
835  slow_energy = slow_gradient = 0;
836  }
837  break;
838  case SPLIT_MARTINI: {
839  // in Martini, the Coulomb switching distance is zero
840  const BigReal COUL_SWITCH = 0.;
841  // Gromacs shifting function
842  const BigReal p1 = 1.;
843  BigReal A1 = p1 * ((p1+1)*COUL_SWITCH-(p1+4)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,2));
844  BigReal B1 = -p1 * ((p1+1)*COUL_SWITCH-(p1+3)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,3));
845  BigReal X1 = 1.0/pow(cutoff,p1)-A1/3.0*pow(cutoff-COUL_SWITCH,3)-B1/4.0*pow(cutoff-COUL_SWITCH,4);
846  BigReal r12 = (r-COUL_SWITCH)*(r-COUL_SWITCH);
847  BigReal r13 = (r-COUL_SWITCH)*(r-COUL_SWITCH)*(r-COUL_SWITCH);
848  BigReal shiftVal = -(A1/3.0)*r13 - (B1/4.0)*r12*r12 - X1;
849  BigReal dShiftVal = -A1*r12 - B1*r13;
850  fast_energy = (1/r) + shiftVal;
851  fast_gradient = -1/(r2) + dShiftVal;
852  scor_energy = scor_gradient = 0;
853  slow_energy = slow_gradient = 0;
854  }
855  break;
856  case SPLIT_C1:
857  // calculate actual energy and gradient
858  slow_energy = 0.5/cutoff * (3.0 - (r2/cutoff2));
859  slow_gradient = -1.0/cutoff2 * (r/cutoff);
860  // calculate scor from slow and corr
861  scor_energy = slow_energy + (corr_energy - 1.0)/r;
862  scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
863  // calculate fast from slow
864  fast_energy = 1.0/r - slow_energy;
865  fast_gradient = -1.0/r2 - slow_gradient;
866  break;
867  case SPLIT_C2:
868  //
869  // Quintic splitting function contributed by
870  // Bruce Berne, Ruhong Zhou, and Joe Morrone
871  //
872  // calculate actual energy and gradient
873  slow_energy = r2/(cutoff*cutoff2) * (6.0 * (r2/cutoff2)
874  - 15.0*(r/cutoff) + 10.0);
875  slow_gradient = r/(cutoff*cutoff2) * (24.0 * (r2/cutoff2)
876  - 45.0 *(r/cutoff) + 20.0);
877  // calculate scor from slow and corr
878  scor_energy = slow_energy + (corr_energy - 1.0)/r;
879  scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
880  // calculate fast from slow
881  fast_energy = 1.0/r - slow_energy;
882  fast_gradient = -1.0/r2 - slow_gradient;
883  break;
884  }
885 
886  // foo_gradient is calculated as ( d foo_energy / d r )
887  // and now divided by 2r to get ( d foo_energy / d r2 )
888 
889  fast_gradient *= 0.5 * r_1;
890  scor_gradient *= 0.5 * r_1;
891  slow_gradient *= 0.5 * r_1;
892 
893  // let modf be 1 if excluded, 1-scale14 if modified, 0 otherwise,
894  // add scor_ - modf * slow_ to slow terms and
895  // add fast_ - modf * fast_ to fast terms.
896 
897  BigReal vdwa_energy, vdwa_gradient;
898  BigReal vdwb_energy, vdwb_gradient;
899 
900  const BigReal r_6 = r_2*r_2*r_2;
901  const BigReal r_12 = r_6*r_6;
902 
903  // Lennard-Jones switching function
904  if ( simParams->vdwForceSwitching ) { // switch force
906 
907  // from Steinbach & Brooks, JCC 15, pgs 667-683, 1994, eqns 10-13
908  if ( r2 > switchOn2 ) {
909  BigReal tmpa = r_6 - cutoff_6;
910  vdwa_energy = k_vdwa * tmpa * tmpa;
911  BigReal tmpb = r_1 * r_2 - cutoff_3;
912  vdwb_energy = k_vdwb * tmpb * tmpb;
913  vdwa_gradient = -6.0 * k_vdwa * tmpa * r_2 * r_6;
914  vdwb_gradient = -3.0 * k_vdwb * tmpb * r_2 * r_2 * r_1;
915  } else {
916  vdwa_energy = r_12 + v_vdwa;
917  vdwb_energy = r_6 + v_vdwb;
918  vdwa_gradient = -6.0 * r_2 * r_12;
919  vdwb_gradient = -3.0 * r_2 * r_6;
920  }
921  } else if ( simParams->martiniSwitching ) { // switching fxn for Martini RBCG
923 
924  BigReal r12 = (r-switchOn)*(r-switchOn); BigReal r13 = (r-switchOn)*(r-switchOn)*(r-switchOn);
925 
926  BigReal p6 = 6;
927  BigReal A6 = p6 * ((p6+1)*switchOn-(p6+4)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,2));
928  BigReal B6 = -p6 * ((p6+1)*switchOn-(p6+3)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,3));
929  BigReal C6 = 1.0/pow(cutoff,p6)-A6/3.0*pow(cutoff-switchOn,3)-B6/4.0*pow(cutoff-switchOn,4);
930 
931  BigReal p12 = 12;
932  BigReal A12 = p12 * ((p12+1)*switchOn-(p12+4)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,2));
933  BigReal B12 = -p12 * ((p12+1)*switchOn-(p12+3)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,3));
934  BigReal C12 = 1.0/pow(cutoff,p12)-A12/3.0*pow(cutoff-switchOn,3)-B12/4.0*pow(cutoff-switchOn,4);
935 
936  BigReal LJshifttempA = -(A12/3)*r13 - (B12/4)*r12*r12 - C12;
937  BigReal LJshifttempB = -(A6/3)*r13 - (B6/4)*r12*r12 - C6;
938  const BigReal shiftValA = // used for Lennard-Jones
939  ( r2 > switchOn2 ? LJshifttempA : -C12);
940  const BigReal shiftValB = // used for Lennard-Jones
941  ( r2 > switchOn2 ? LJshifttempB : -C6);
942 
943  BigReal LJdshifttempA = -A12*r12 - B12*r13;
944  BigReal LJdshifttempB = -A6*r12 - B6*r13;
945  const BigReal dshiftValA = // used for Lennard-Jones
946  ( r2 > switchOn2 ? LJdshifttempA*0.5*r_1 : 0 );
947  const BigReal dshiftValB = // used for Lennard-Jones
948  ( r2 > switchOn2 ? LJdshifttempB*0.5*r_1 : 0 );
949 
950 
951 
952 
953  //have not addressed r > cutoff
954 
955  // dshiftValA*= 0.5*r_1;
956  // dshiftValB*= 0.5*r_1;
957 
958  vdwa_energy = r_12 + shiftValA;
959  vdwb_energy = r_6 + shiftValB;
960 
961  vdwa_gradient = -6/pow(r,14) + dshiftValA ;
962  vdwb_gradient = -3/pow(r,8) + dshiftValB;
963 
964  } else { // switch energy
966 
967  const BigReal c2 = cutoff2-r2;
968  const BigReal c4 = c2*(c3-2.0*c2);
969  const BigReal switchVal = // used for Lennard-Jones
970  ( r2 > switchOn2 ? c2*c4*c1 : 1.0 );
971  const BigReal dSwitchVal = // d switchVal / d r2
972  ( r2 > switchOn2 ? 2*c1*(c2*c2-c4) : 0.0 );
973 
974  vdwa_energy = switchVal * r_12;
975  vdwb_energy = switchVal * r_6;
976 
977  vdwa_gradient = ( dSwitchVal - 6.0 * switchVal * r_2 ) * r_12;
978  vdwb_gradient = ( dSwitchVal - 3.0 * switchVal * r_2 ) * r_6;
979  }
980 
981 
982 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
983  if ( knl_table ) {
984  knl_fast_ener_table[i] = -1.*fast_energy;
985  knl_fast_grad_table[i] = -2.*fast_gradient;
986  knl_scor_ener_table[i] = -1.*scor_energy;
987  knl_scor_grad_table[i] = -2.*scor_gradient;
988  knl_slow_ener_table[i] = (-1.*scor_energy - (scale14-1)*slow_energy) /
989  scale14;
990  knl_slow_grad_table[i] = (-2.*scor_gradient -
991  (scale14-1)*2.*slow_gradient) / scale14;
992  knl_excl_ener_table[i] = slow_energy - scor_energy;
993  knl_excl_grad_table[i] = 2.*(slow_gradient - scor_gradient);
994 #ifdef NAMD_KNL
995  knl_corr_ener_table[i] = -1.*(fast_energy + scor_energy);
996  knl_corr_grad_table[i] = -2.*(fast_gradient + scor_gradient);
997 #endif
998  if ( i == nn-1 ) {
999  knl_fast_ener_table[nn] = knl_fast_ener_table[i];
1000  knl_fast_grad_table[nn] = knl_fast_grad_table[i];
1001  knl_scor_ener_table[nn] = knl_scor_ener_table[i];
1002  knl_scor_grad_table[nn] = knl_scor_grad_table[i];
1003  knl_slow_ener_table[nn] = knl_slow_ener_table[i];
1004  knl_slow_grad_table[nn] = knl_slow_grad_table[i];
1005  knl_excl_ener_table[nn] = knl_excl_ener_table[i];
1006  knl_excl_grad_table[nn] = knl_excl_grad_table[i];
1007 #ifdef NAMD_KNL
1008  knl_corr_ener_table[nn] = knl_corr_ener_table[i];
1009  knl_corr_grad_table[nn] = knl_corr_grad_table[i];
1010 #endif
1011  }
1012  } else {
1013 #endif
1014  *(fast_i++) = fast_energy;
1015  *(fast_i++) = fast_gradient;
1016  *(fast_i++) = 0;
1017  *(fast_i++) = 0;
1018  *(scor_i++) = scor_energy;
1019  *(scor_i++) = scor_gradient;
1020  *(scor_i++) = 0;
1021  *(scor_i++) = 0;
1022  *(slow_i++) = slow_energy;
1023  *(slow_i++) = slow_gradient;
1024  *(slow_i++) = 0;
1025  *(slow_i++) = 0;
1026  *(vdwa_i++) = vdwa_energy;
1027  *(vdwa_i++) = vdwa_gradient;
1028  *(vdwa_i++) = 0;
1029  *(vdwa_i++) = 0;
1030  *(vdwb_i++) = vdwb_energy;
1031  *(vdwb_i++) = vdwb_gradient;
1032  *(vdwb_i++) = 0;
1033  *(vdwb_i++) = 0;
1034  *(r2_i++) = r2 + r2_delta;
1035 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
1036  }
1037 #endif
1038 
1039  }
1040 #if defined(NAMD_KNL) || defined(NAMD_AVXTILES)
1041  } // knl_table loop
1042 #endif
1043 
1044  if ( ! r2_delta_i ) {
1045  NAMD_bug("Failed to find table entry for r2 == r2_limit\n");
1046  }
1047  if ( r2_table[r2_delta_i] > r2_limit + r2_delta ) {
1048  NAMD_bug("Found bad table entry for r2 == r2_limit\n");
1049  }
1050 
1051  int j;
1052  const char *table_name = "XXXX";
1053  int smooth_short = 0;
1054  for ( j=0; j<5; ++j ) {
1055  BigReal *t0 = 0;
1056  switch (j) {
1057  case 0:
1058  t0 = fast_table;
1059  table_name = "FAST";
1060  smooth_short = 1;
1061  break;
1062  case 1:
1063  t0 = scor_table;
1064  table_name = "SCOR";
1065  smooth_short = 0;
1066  break;
1067  case 2:
1068  t0 = slow_table;
1069  table_name = "SLOW";
1070  smooth_short = 0;
1071  break;
1072  case 3:
1073  t0 = vdwa_table;
1074  table_name = "VDWA";
1075  smooth_short = 1;
1076  break;
1077  case 4:
1078  t0 = vdwb_table;
1079  table_name = "VDWB";
1080  smooth_short = 1;
1081  break;
1082  }
1083  // patch up data for i=0
1084  t0[0] = t0[4] - t0[5] * ( r2_delta / 64.0 ); // energy
1085  t0[1] = t0[5]; // gradient
1086  t0[2] = 0;
1087  t0[3] = 0;
1088  if ( smooth_short ) {
1089  BigReal energy0 = t0[4*r2_delta_i];
1090  BigReal gradient0 = t0[4*r2_delta_i+1];
1091  BigReal r20 = r2_table[r2_delta_i];
1092  t0[0] = energy0 - gradient0 * (r20 - r2_table[0]); // energy
1093  t0[1] = gradient0; // gradient
1094  }
1095  BigReal *t;
1096  for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
1097  BigReal x = ( r2_delta * ( 1 << (i/64) ) ) / 64.0;
1098  if ( r2_table[i+1] != r2_table[i] + x ) {
1099  NAMD_bug("Bad table delta calculation.\n");
1100  }
1101  if ( smooth_short && i+1 < r2_delta_i ) {
1102  BigReal energy0 = t0[4*r2_delta_i];
1103  BigReal gradient0 = t0[4*r2_delta_i+1];
1104  BigReal r20 = r2_table[r2_delta_i];
1105  t[4] = energy0 - gradient0 * (r20 - r2_table[i+1]); // energy
1106  t[5] = gradient0; // gradient
1107  }
1108  BigReal v1 = t[0];
1109  BigReal g1 = t[1];
1110  BigReal v2 = t[4];
1111  BigReal g2 = t[5];
1112  // explicit formulas for v1 + g1 x + c x^2 + d x^3
1113  BigReal c = ( 3.0 * (v2 - v1) - x * (2.0 * g1 + g2) ) / ( x * x );
1114  BigReal d = ( -2.0 * (v2 - v1) + x * (g1 + g2) ) / ( x * x * x );
1115  // since v2 - v1 is imprecise, we refine c and d numerically
1116  // important because we need accurate forces (more than energies!)
1117  for ( int k=0; k < 2; ++k ) {
1118  BigReal dv = (v1 - v2) + ( ( d * x + c ) * x + g1 ) * x;
1119  BigReal dg = (g1 - g2) + ( 3.0 * d * x + 2.0 * c ) * x;
1120  c -= ( 3.0 * dv - x * dg ) / ( x * x );
1121  d -= ( -2.0 * dv + x * dg ) / ( x * x * x );
1122  }
1123  // store in the array;
1124  t[2] = c; t[3] = d;
1125  }
1126 
1127  if ( ! CkMyPe() ) {
1128  BigReal dvmax = 0;
1129  BigReal dgmax = 0;
1130  BigReal dvmax_r = 0;
1131  BigReal dgmax_r = 0;
1132  BigReal fdvmax = 0;
1133  BigReal fdgmax = 0;
1134  BigReal fdvmax_r = 0;
1135  BigReal fdgmax_r = 0;
1136  BigReal dgcdamax = 0;
1137  BigReal dgcdimax = 0;
1138  BigReal dgcaimax = 0;
1139  BigReal dgcdamax_r = 0;
1140  BigReal dgcdimax_r = 0;
1141  BigReal dgcaimax_r = 0;
1142  BigReal fdgcdamax = 0;
1143  BigReal fdgcdimax = 0;
1144  BigReal fdgcaimax = 0;
1145  BigReal fdgcdamax_r = 0;
1146  BigReal fdgcdimax_r = 0;
1147  BigReal fdgcaimax_r = 0;
1148  BigReal gcm = fabs(t0[1]); // gradient magnitude running average
1149  for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
1150  const BigReal r2_base = r2_delta * ( 1 << (i/64) );
1151  const BigReal r2_del = r2_base / 64.0;
1152  const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
1153  const BigReal r = sqrt(r2);
1154  if ( r > cutoff ) break;
1155  BigReal x = r2_del;
1156  BigReal dv = ( ( t[3] * x + t[2] ) * x + t[1] ) * x + t[0] - t[4];
1157  BigReal dg = ( 3.0 * t[3] * x + 2.0 * t[2] ) * x + t[1] - t[5];
1158  if ( t[4] != 0. && fabs(dv/t[4]) > fdvmax ) {
1159  fdvmax = fabs(dv/t[4]); fdvmax_r = r;
1160  }
1161  if ( fabs(dv) > dvmax ) {
1162  dvmax = fabs(dv); dvmax_r = r;
1163  }
1164  if ( t[5] != 0. && fabs(dg/t[5]) > fdgmax ) {
1165  fdgmax = fabs(dg/t[5]); fdgmax_r = r;
1166  }
1167  if ( fabs(dg) > dgmax ) {
1168  dgmax = fabs(dg); dgmax_r = r;
1169  }
1170  BigReal gcd = (t[4] - t[0]) / x; // centered difference gradient
1171  BigReal gcd_prec = (fabs(t[0]) + fabs(t[4])) * 1.e-15 / x; // roundoff
1172  gcm = 0.9 * gcm + 0.1 * fabs(t[5]); // magnitude running average
1173  BigReal gca = 0.5 * (t[1] + t[5]); // centered average gradient
1174  BigReal gci = ( 0.75 * t[3] * x + t[2] ) * x + t[1]; // interpolated
1175  BigReal rc = sqrt(r2 + 0.5 * x);
1176  BigReal dgcda = gcd - gca;
1177  if ( dgcda != 0. && fabs(dgcda) < gcd_prec ) {
1178  // CkPrintf("ERROR %g < PREC %g AT %g AVG VAL %g\n", dgcda, gcd_prec, rc, gca);
1179  dgcda = 0.;
1180  }
1181  BigReal dgcdi = gcd - gci;
1182  if ( dgcdi != 0. && fabs(dgcdi) < gcd_prec ) {
1183  // CkPrintf("ERROR %g < PREC %g AT %g INT VAL %g\n", dgcdi, gcd_prec, rc, gci);
1184  dgcdi = 0.;
1185  }
1186  BigReal dgcai = gca - gci;
1187  if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcda/gcm) > fdgcdamax ) {
1188  fdgcdamax = fabs(dgcda/gcm); fdgcdamax_r = rc;
1189  }
1190  if ( fabs(dgcda) > fdgcdamax ) {
1191  dgcdamax = fabs(dgcda); dgcdamax_r = rc;
1192  }
1193  if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcdi/gcm) > fdgcdimax ) {
1194  fdgcdimax = fabs(dgcdi/gcm); fdgcdimax_r = rc;
1195  }
1196  if ( fabs(dgcdi) > fdgcdimax ) {
1197  dgcdimax = fabs(dgcdi); dgcdimax_r = rc;
1198  }
1199  if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcai/gcm) > fdgcaimax ) {
1200  fdgcaimax = fabs(dgcai/gcm); fdgcaimax_r = rc;
1201  }
1202  if ( fabs(dgcai) > fdgcaimax ) {
1203  dgcaimax = fabs(dgcai); dgcaimax_r = rc;
1204  }
1205 #if 0
1206  CkPrintf("TABLE %s %g %g %g %g\n",table_name,rc,dgcda/gcm,dgcda,gci);
1207  if (dv != 0.) CkPrintf("TABLE %d ENERGY ERROR %g AT %g (%d)\n",j,dv,r,i);
1208  if (dg != 0.) CkPrintf("TABLE %d FORCE ERROR %g AT %g (%d)\n",j,dg,r,i);
1209 #endif
1210  }
1211  if ( dvmax != 0.0 ) {
1212  iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
1213  " TABLE ENERGY: " << dvmax << " AT " << dvmax_r << "\n" << endi;
1214  }
1215  if ( fdvmax != 0.0 ) {
1216  iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
1217  " TABLE ENERGY: " << fdvmax << " AT " << fdvmax_r << "\n" << endi;
1218  }
1219  if ( dgmax != 0.0 ) {
1220  iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
1221  " TABLE FORCE: " << dgmax << " AT " << dgmax_r << "\n" << endi;
1222  }
1223  if ( fdgmax != 0.0 ) {
1224  iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
1225  " TABLE FORCE: " << fdgmax << " AT " << fdgmax_r << "\n" << endi;
1226  }
1227  if (fdgcdamax != 0.0 ) {
1228  iout << iINFO << "INCONSISTENCY IN " << table_name <<
1229  " TABLE ENERGY VS FORCE: " << fdgcdamax << " AT " << fdgcdamax_r << "\n" << endi;
1230  if ( fdgcdamax > 0.1 ) {
1231  iout << iERROR << "\n";
1232  iout << iERROR << "CALCULATED " << table_name <<
1233  " FORCE MAY NOT MATCH ENERGY! POSSIBLE BUG!\n";
1234  iout << iERROR << "\n";
1235  }
1236  }
1237  if (0 && fdgcdimax != 0.0 ) {
1238  iout << iINFO << "INCONSISTENCY IN " << table_name <<
1239  " TABLE ENERGY VS FORCE: " << fdgcdimax << " AT " << fdgcdimax_r << "\n" << endi;
1240  }
1241  if ( 0 && fdgcaimax != 0.0 ) {
1242  iout << iINFO << "INCONSISTENCY IN " << table_name <<
1243  " TABLE AVG VS INT FORCE: " << fdgcaimax << " AT " << fdgcaimax_r << "\n" << endi;
1244  }
1245  }
1246 
1247  }
1248 
1249  for ( i=0; i<4*n; ++i ) {
1250  corr_table[i] = fast_table[i] + scor_table[i];
1251  full_table[i] = fast_table[i] + slow_table[i];
1252  }
1253 
1254 #if 0
1255  for ( i=0; i<n; ++i ) {
1256  for ( int j=0; j<4; ++j ) {
1257  table_short[16*i+6-2*j] = table_noshort[16*i+6-2*j] = vdwa_table[4*i+j];
1258  table_short[16*i+7-2*j] = table_noshort[16*i+7-2*j] = vdwb_table[4*i+j];
1259  table_short[16*i+8+3-j] = fast_table[4*i+j];
1260  table_short[16*i+12+3-j] = scor_table[4*i+j];
1261  table_noshort[16*i+8+3-j] = corr_table[4*i+j];
1262  table_noshort[16*i+12+3-j] = full_table[4*i+j];
1263  }
1264  }
1265 #endif
1266 
1267  for ( i=0; i<n; ++i ) {
1268  table_short[16*i+0] = table_noshort[16*i+0] = -6.*vdwa_table[4*i+3];
1269  table_short[16*i+1] = table_noshort[16*i+1] = -4.*vdwa_table[4*i+2];
1270  table_short[16*i+2] = table_noshort[16*i+2] = -2.*vdwa_table[4*i+1];
1271  table_short[16*i+3] = table_noshort[16*i+3] = -1.*vdwa_table[4*i+0];
1272 
1273  table_short[16*i+4] = table_noshort[16*i+4] = -6.*vdwb_table[4*i+3];
1274  table_short[16*i+5] = table_noshort[16*i+5] = -4.*vdwb_table[4*i+2];
1275  table_short[16*i+6] = table_noshort[16*i+6] = -2.*vdwb_table[4*i+1];
1276  table_short[16*i+7] = table_noshort[16*i+7] = -1.*vdwb_table[4*i+0];
1277 
1278  table_short[16*i+8] = -6.*fast_table[4*i+3];
1279  table_short[16*i+9] = -4.*fast_table[4*i+2];
1280  table_short[16*i+10] = -2.*fast_table[4*i+1];
1281  table_short[16*i+11] = -1.*fast_table[4*i+0];
1282 
1283  table_noshort[16*i+8] = -6.*corr_table[4*i+3];
1284  table_noshort[16*i+9] = -4.*corr_table[4*i+2];
1285  table_noshort[16*i+10] = -2.*corr_table[4*i+1];
1286  table_noshort[16*i+11] = -1.*corr_table[4*i+0];
1287 
1288  table_short[16*i+12] = -6.*scor_table[4*i+3];
1289  table_short[16*i+13] = -4.*scor_table[4*i+2];
1290  table_short[16*i+14] = -2.*scor_table[4*i+1];
1291  table_short[16*i+15] = -1.*scor_table[4*i+0];
1292 
1293  table_noshort[16*i+12] = -6.*full_table[4*i+3];
1294  table_noshort[16*i+13] = -4.*full_table[4*i+2];
1295  table_noshort[16*i+14] = -2.*full_table[4*i+1];
1296  table_noshort[16*i+15] = -1.*full_table[4*i+0];
1297  }
1298 
1299 #if 0
1300  char fname[100];
1301  sprintf(fname,"/tmp/namd.table.pe%d.dat",CkMyPe());
1302  FILE *f = fopen(fname,"w");
1303  for ( i=0; i<(n-1); ++i ) {
1304  const BigReal r2_base = r2_delta * ( 1 << (i/64) );
1305  const BigReal r2_del = r2_base / 64.0;
1306  const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
1307  BigReal *t;
1308  if ( r2 + r2_delta != r2_table[i] ) fprintf(f,"r2 error! ");
1309  fprintf(f,"%g",r2);
1310  t = fast_table + 4*i;
1311  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1312  t = scor_table + 4*i;
1313  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1314  t = slow_table + 4*i;
1315  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1316  t = corr_table + 4*i;
1317  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1318  t = full_table + 4*i;
1319  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1320  t = vdwa_table + 4*i;
1321  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1322  t = vdwb_table + 4*i;
1323  fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
1324  fprintf(f,"\n");
1325  }
1326  fclose(f);
1327 #endif
1328 
1329  //Flip slow table to match table_four_i
1330  for ( i=0; i<n; ++i ) {
1331  BigReal tmp0, tmp1, tmp2, tmp3;
1332  tmp0 = slow_table [i*4 + 0];
1333  tmp1 = slow_table [i*4 + 1];
1334  tmp2 = slow_table [i*4 + 2];
1335  tmp3 = slow_table [i*4 + 3];
1336 
1337  slow_table [i*4 + 0] = tmp3;
1338  slow_table [i*4 + 1] = tmp2;
1339  slow_table [i*4 + 2] = tmp1;
1340  slow_table [i*4 + 3] = tmp0;
1341  }
1342 
1343 #ifdef NAMD_CUDA
1344  if (!simParams->useCUDA2) {
1346  }
1347 #endif
1348 
1349  #ifdef NAMD_MIC
1350  send_build_mic_force_table();
1351  #endif
1352 }
1353 
static Node * Object()
Definition: Node.h:86
static void calc_pair_fullelect_pprof(nonbonded *)
static BigReal * fast_table
#define SCALED14
Definition: SimParameters.h:43
static void calc_pair_energy_merge_fullelect_tabener(nonbonded *)
static void calc_self_energy_go(nonbonded *)
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
static void calc_self_energy_merge_fullelect_ti(nonbonded *)
static BigReal * scor_table
static void calc_pair_les(nonbonded *)
static void calc_self_energy_slow_fullelect_fep(nonbonded *)
static void calc_pair_merge_fullelect(nonbonded *)
int columnsize
Definition: Parameters.h:253
static void calc_self_energy_merge_fullelect_int(nonbonded *)
static void calc_self_energy_slow_fullelect(nonbonded *)
#define XPLOR
Definition: SimParameters.h:55
static void calc_self_energy_merge_fullelect_tabener(nonbonded *)
Bool vdwGeometricSigma
static void calc_self_energy_slow_fullelect_ti(nonbonded *)
static void calc_self_energy_fullelect_int(nonbonded *)
static void submitReductionData(BigReal *, SubmitReduction *)
static void calc_self_ti(nonbonded *)
static void calc_self_fullelect_pprof(nonbonded *)
static void calc_self_tabener(nonbonded *)
static void calc_self_energy_fullelect(nonbonded *)
static void calc_pair_slow_fullelect_ti(nonbonded *)
#define SHARP
Definition: SimParameters.h:54
static void(* calcSelf)(nonbonded *)
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:32
static __global__ void(const patch_pair *patch_pairs, const atom *atoms, const atom_param *atom_params, const int *vdw_types, unsigned int *plist, float4 *tmpforces, float4 *slow_tmpforces, float4 *forces, float4 *slow_forces, float *tmpvirials, float *slow_tmpvirials, float *virials, float *slow_virials, unsigned int *global_counters, int *force_ready_queue, const unsigned int *overflow_exclusions, const int npatches, const int block_begin, const int total_block_count, int *block_order, exclmask *exclmasks, const int lj_table_size, const float3 lata, const float3 latb, const float3 latc, const float cutoff2, const float plcutoff2, const int doSlow)
static void calc_pair_merge_fullelect_tabener(nonbonded *)
static void calc_self_les(nonbonded *)
const BigReal A
SimParameters * simParameters
Definition: Node.h:178
BigReal nonbondedScaling
static void calc_pair_merge_fullelect_go(nonbonded *)
#define SPLIT_NONE
static void calc_pair_energy_merge_fullelect(nonbonded *)
static void calc_pair_slow_fullelect_pprof(nonbonded *)
static const Molecule * mol
static void calc_pair_energy_fullelect_fep(nonbonded *)
static void calc_self_energy_fullelect_pprof(nonbonded *)
float Real
Definition: common.h:109
BigReal & item(int i)
Definition: ReductionMgr.h:312
static void calc_self_energy_fullelect_les(nonbonded *)
static void calc_self_energy_fullelect_ti(nonbonded *)
BigReal limitDist
static void calc_self_fullelect(nonbonded *)
#define VDW_SWITCH_MODE_MARTINI
static void calc_self_energy_ti(nonbonded *)
static BigReal * vdwa_table
static void calc_self_energy_int(nonbonded *)
static void submitPressureProfileData(BigReal *, SubmitReduction *)
static BigReal pressureProfileThickness
int get_table_dim() const
Definition: LJTable.h:44
static void calc_pair_slow_fullelect(nonbonded *)
static void calc_pair(nonbonded *)
static void calc_self_go(nonbonded *)
#define iout
Definition: InfoStream.h:87
static void(* calcMergePair)(nonbonded *)
int pressureProfileSlabs
static void calc_pair_energy_tabener(nonbonded *)
static void calc_self_slow_fullelect_pprof(nonbonded *)
static void calc_pair_energy_slow_fullelect_go(nonbonded *)
static void calc_self_energy_merge_fullelect_fep(nonbonded *)
static void(* calcMergePairEnergy)(nonbonded *)
static void calc_pair_energy_go(nonbonded *)
static void calc_pair_merge_fullelect_ti(nonbonded *)
#define C1
Definition: SimParameters.h:56
Bool pairInteractionOn
static BigReal * full_table
static void calc_pair_merge_fullelect_les(nonbonded *)
static void calc_pair_slow_fullelect_go(nonbonded *)
#define SPLIT_MARTINI
static void calc_pair_energy_merge_fullelect_fep(nonbonded *)
BigReal switchingDist
BigReal drudeNbtholeCut
static void(* calcSlowPairEnergy)(nonbonded *)
void send_build_cuda_force_table()
static void(* calcPair)(nonbonded *)
void add(int nitems, const BigReal *arr)
Definition: ReductionMgr.h:320
static void(* calcSlowPair)(nonbonded *)
#define SPLIT_SHIFT
static void calc_pair_energy_fullelect_go(nonbonded *)
static void calc_pair_energy_slow_fullelect_tabener(nonbonded *)
static void calc_self_slow_fullelect(nonbonded *)
static void calc_self_energy_les(nonbonded *)
static void calc_pair_energy_fep(nonbonded *)
static void(* calcMergeSelfEnergy)(nonbonded *)
static void calc_pair_tabener(nonbonded *)
static void calc_pair_energy_pprof(nonbonded *)
static BigReal * table_noshort
static void calc_self_energy_slow_fullelect_tabener(nonbonded *)
void NAMD_bug(const char *err_msg)
Definition: common.C:123
static void calc_self_energy_fullelect_tabener(nonbonded *)
static void calc_pair_energy_int(nonbonded *)
static BigReal pressureProfileMin
static void calc_self_slow_fullelect_tabener(nonbonded *)
static void calc_pair_energy_slow_fullelect_les(nonbonded *)
static void calc_pair_energy_ti(nonbonded *)
static void calc_pair_energy_slow_fullelect_fep(nonbonded *)
static void calc_self_slow_fullelect_ti(nonbonded *)
static void calc_pair_slow_fullelect_les(nonbonded *)
GoChoices goMethod
static void calc_pair_fullelect_ti(nonbonded *)
static void calc_pair_energy(nonbonded *)
static BigReal * table_ener
#define C2
Definition: SimParameters.h:57
int Bool
Definition: common.h:133
static void calc_pair_pprof(nonbonded *)
#define ADD_VECTOR(R, RL, D, DL)
Definition: ReductionMgr.h:22
static void calc_pair_energy_fullelect_ti(nonbonded *)
static void calc_self_energy_merge_fullelect_les(nonbonded *)
static void calc_self_slow_fullelect_go(nonbonded *)
static void calc_pair_energy_merge_fullelect_go(nonbonded *)
static void calc_error(nonbonded *)
static void calc_self_energy_fullelect_go(nonbonded *)
BigReal PMEEwaldCoefficient
int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *)
Definition: Parameters.C:3680
static void calc_self_fullelect_tabener(nonbonded *)
static void(* calcFullSelf)(nonbonded *)
void NAMD_die(const char *err_msg)
Definition: common.C:83
static void calc_self_merge_fullelect_ti(nonbonded *)
static void(* calcSlowSelf)(nonbonded *)
#define SPLIT_C1
static void calc_self_fullelect_les(nonbonded *)
static void calc_self_merge_fullelect_go(nonbonded *)
int get_num_vdw_params(void)
Definition: Parameters.h:531
static void calc_self_energy_merge_fullelect_go(nonbonded *)
static void calc_self_energy_merge_fullelect(nonbonded *)
Bool vdwForceSwitching
static void select(void)
static void calc_pair_energy_fullelect_int(nonbonded *)
static void calc_pair_energy_merge_fullelect_pprof(nonbonded *)
static BigReal * lambda_table
static void calc_self_energy_fep(nonbonded *)
static void(* calcSlowSelfEnergy)(nonbonded *)
static void calc_pair_energy_slow_fullelect_pprof(nonbonded *)
Parameters * parameters
Definition: Node.h:177
static BigReal * slow_table
static void(* calcSelfEnergy)(nonbonded *)
static void(* calcPairEnergy)(nonbonded *)
static void calc_pair_fullelect_les(nonbonded *)
static void calc_pair_energy_les(nonbonded *)
int pressureProfileAtomTypes
#define simParams
Definition: Output.C:127
static void calc_self_merge_fullelect(nonbonded *)
static void calc_self_energy_tabener(nonbonded *)
static void calc_pair_energy_merge_fullelect_int(nonbonded *)
static void calc_self_energy_slow_fullelect_pprof(nonbonded *)
static void calc_self_energy_pprof(nonbonded *)
static void calc_pair_energy_fullelect_pprof(nonbonded *)
static void calc_pair_fullelect_go(nonbonded *)
static void calc_pair_energy_slow_fullelect(nonbonded *)
static void calc_pair_energy_fullelect_les(nonbonded *)
static BigReal * vdwb_table
static void calc_pair_energy_fullelect_tabener(nonbonded *)
#define SPLIT_C2
BigReal * table_ener
Definition: Parameters.h:254
static void(* calcFullPair)(nonbonded *)
#define VDW_SWITCH_MODE_FORCE
static void calc_pair_merge_fullelect_pprof(nonbonded *)
static const LJTable * ljTable
static void calc_self_energy_merge_fullelect_pprof(nonbonded *)
const BigReal B
#define VDW_SWITCH_MODE_ENERGY
Bool pressureProfileOn
static void(* calcFullPairEnergy)(nonbonded *)
static BigReal * corr_table
static BigReal * table_short
Bool pairInteractionSelf
static void calc_self(nonbonded *)
#define SPOLY(pg, pdg, ra, split)
Definition: msm_defn.h:140
static void calc_self_fullelect_ti(nonbonded *)
BigReal dielectric
static void calc_pair_fullelect_tabener(nonbonded *)
static void calc_pair_energy_merge_fullelect_ti(nonbonded *)
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
static BigReal * table_alloc
static void calc_self_energy_fullelect_fep(nonbonded *)
Bool tabulatedEnergies
static void calc_self_merge_fullelect_tabener(nonbonded *)
static BigReal alchVdwShiftCoeff
static void calc_pair_go(nonbonded *)
gridSize x
static void calc_pair_ti(nonbonded *)
ExclusionSettings exclude
static void calc_self_energy(nonbonded *)
void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, Real *epsilon14, Index index)
Definition: Parameters.h:497
static void calc_self_fullelect_go(nonbonded *)
static void calc_pair_energy_merge_fullelect_les(nonbonded *)
static void calc_self_pprof(nonbonded *)
static void calc_self_merge_fullelect_les(nonbonded *)
Molecule * molecule
Definition: Node.h:176
static void calc_pair_slow_fullelect_tabener(nonbonded *)
static void calc_pair_energy_fullelect(nonbonded *)
static void calc_self_merge_fullelect_pprof(nonbonded *)
static void calc_self_energy_slow_fullelect_go(nonbonded *)
static void calc_pair_energy_slow_fullelect_ti(nonbonded *)
static void(* calcMergeSelf)(nonbonded *)
static void calc_self_energy_slow_fullelect_les(nonbonded *)
double BigReal
Definition: common.h:114
static void(* calcFullSelfEnergy)(nonbonded *)
BigReal alchVdwShiftCoeff
static void calc_pair_fullelect(nonbonded *)
static void calc_self_slow_fullelect_les(nonbonded *)