ComputeNonbondedUtil.C

Go to the documentation of this file.
00001 
00007 /*
00008    Common operations for ComputeNonbonded classes
00009 */
00010 
00011 #ifdef WIN32
00012 extern "C" {
00013   double erfc(double);
00014 }
00015 #endif
00016 
00017 #include "InfoStream.h"
00018 #include "ComputeNonbondedUtil.h"
00019 #include "SimParameters.h"
00020 #include "Node.h"
00021 #include "Molecule.h"
00022 #include "LJTable.h"
00023 #include "ReductionMgr.h"
00024 #include "Parameters.h"
00025 #include "MsmMacros.h"
00026 #include <stdio.h>
00027 
00028 #ifdef NAMD_CUDA
00029   void send_build_cuda_force_table();
00030 #endif
00031 
00032 #ifdef NAMD_MIC
00033   extern void send_build_mic_force_table();
00034 #endif
00035 
00036 Bool            ComputeNonbondedUtil::commOnly;
00037 Bool            ComputeNonbondedUtil::fixedAtomsOn;
00038 Bool            ComputeNonbondedUtil::qmForcesOn;
00039 BigReal         ComputeNonbondedUtil::cutoff;
00040 BigReal         ComputeNonbondedUtil::cutoff2;
00041 float           ComputeNonbondedUtil::cutoff2_f;
00042 BigReal         ComputeNonbondedUtil::dielectric_1;
00043 const LJTable*  ComputeNonbondedUtil::ljTable = 0;
00044 const Molecule* ComputeNonbondedUtil::mol;
00045 BigReal         ComputeNonbondedUtil::r2_delta;
00046 BigReal         ComputeNonbondedUtil::r2_delta_1;
00047 int             ComputeNonbondedUtil::r2_delta_exp;
00048 int             ComputeNonbondedUtil::rowsize;
00049 int             ComputeNonbondedUtil::columnsize;
00050 BigReal*        ComputeNonbondedUtil::table_alloc = 0;
00051 BigReal*        ComputeNonbondedUtil::table_ener = 0;
00052 BigReal*        ComputeNonbondedUtil::table_short;
00053 BigReal*        ComputeNonbondedUtil::table_noshort;
00054 BigReal*        ComputeNonbondedUtil::fast_table;
00055 BigReal*        ComputeNonbondedUtil::scor_table;
00056 BigReal*        ComputeNonbondedUtil::slow_table;
00057 BigReal*        ComputeNonbondedUtil::corr_table;
00058 BigReal*        ComputeNonbondedUtil::full_table;
00059 BigReal*        ComputeNonbondedUtil::vdwa_table;
00060 BigReal*        ComputeNonbondedUtil::vdwb_table;
00061 BigReal*        ComputeNonbondedUtil::r2_table;
00062 int ComputeNonbondedUtil::table_length;
00063 #if defined(NAMD_MIC)
00064   BigReal*      ComputeNonbondedUtil::mic_table_base_ptr;
00065   int           ComputeNonbondedUtil::mic_table_n;
00066   int           ComputeNonbondedUtil::mic_table_n_16;
00067 #endif
00068 #ifdef NAMD_KNL
00069 float*          ComputeNonbondedUtil::knl_table_alloc;
00070 float*          ComputeNonbondedUtil::knl_fast_ener_table;
00071 float*          ComputeNonbondedUtil::knl_fast_grad_table;
00072 float*          ComputeNonbondedUtil::knl_scor_ener_table;
00073 float*          ComputeNonbondedUtil::knl_scor_grad_table;
00074 float*          ComputeNonbondedUtil::knl_slow_ener_table;
00075 float*          ComputeNonbondedUtil::knl_slow_grad_table;
00076 float*          ComputeNonbondedUtil::knl_corr_ener_table;
00077 float*          ComputeNonbondedUtil::knl_corr_grad_table;
00078 float*          ComputeNonbondedUtil::knl_full_ener_table;
00079 float*          ComputeNonbondedUtil::knl_full_grad_table;
00080 #endif
00081 BigReal         ComputeNonbondedUtil::scaling;
00082 BigReal         ComputeNonbondedUtil::scale14;
00083 BigReal         ComputeNonbondedUtil::switchOn;
00084 BigReal         ComputeNonbondedUtil::switchOn_1;
00085 BigReal         ComputeNonbondedUtil::switchOn2;
00086 BigReal         ComputeNonbondedUtil::v_vdwa;
00087 BigReal         ComputeNonbondedUtil::v_vdwb;
00088 BigReal         ComputeNonbondedUtil::k_vdwa;
00089 BigReal         ComputeNonbondedUtil::k_vdwb;
00090 BigReal         ComputeNonbondedUtil::cutoff_3;
00091 BigReal         ComputeNonbondedUtil::cutoff_6;
00092 float           ComputeNonbondedUtil::v_vdwa_f;
00093 float           ComputeNonbondedUtil::v_vdwb_f;
00094 float           ComputeNonbondedUtil::k_vdwa_f;
00095 float           ComputeNonbondedUtil::k_vdwb_f;
00096 float           ComputeNonbondedUtil::cutoff_3_f;
00097 float           ComputeNonbondedUtil::cutoff_6_f;
00098 float           ComputeNonbondedUtil::switchOn_f;
00099 float           ComputeNonbondedUtil::A6_f;
00100 float           ComputeNonbondedUtil::B6_f;
00101 float           ComputeNonbondedUtil::C6_f;
00102 float           ComputeNonbondedUtil::A12_f;
00103 float           ComputeNonbondedUtil::B12_f;
00104 float           ComputeNonbondedUtil::C12_f;
00105 BigReal         ComputeNonbondedUtil::c0;
00106 BigReal         ComputeNonbondedUtil::c1;
00107 BigReal         ComputeNonbondedUtil::c3;
00108 BigReal         ComputeNonbondedUtil::c5;
00109 BigReal         ComputeNonbondedUtil::c6;
00110 BigReal         ComputeNonbondedUtil::c7;
00111 BigReal         ComputeNonbondedUtil::c8;
00112 // BigReal         ComputeNonbondedUtil::d0;
00113 // fepb
00114 Bool      ComputeNonbondedUtil::alchFepOn;
00115 Bool      ComputeNonbondedUtil::alchThermIntOn;
00116 Bool      ComputeNonbondedUtil::Fep_WCA_repuOn;
00117 Bool      ComputeNonbondedUtil::Fep_WCA_dispOn;
00118 Bool      ComputeNonbondedUtil::Fep_ElecOn;
00119 Bool      ComputeNonbondedUtil::Fep_Wham; 
00120 BigReal   ComputeNonbondedUtil::WCA_rcut1;
00121 BigReal   ComputeNonbondedUtil::WCA_rcut2;
00122 BigReal   ComputeNonbondedUtil::WCA_rcut3;
00123 BigReal   ComputeNonbondedUtil::alchLambda2;
00124 BigReal   ComputeNonbondedUtil::alchRepLambda;
00125 BigReal   ComputeNonbondedUtil::alchDispLambda;
00126 BigReal   ComputeNonbondedUtil::alchElecLambda;
00127 BigReal   ComputeNonbondedUtil::alchVdwShiftCoeff;
00128 Bool      ComputeNonbondedUtil::vdwForceSwitching;
00129 Bool      ComputeNonbondedUtil::LJcorrection;
00130 Bool      ComputeNonbondedUtil::alchDecouple;
00131 //fepe
00132 Bool      ComputeNonbondedUtil::lesOn;
00133 int       ComputeNonbondedUtil::lesFactor;
00134 BigReal   ComputeNonbondedUtil::lesScaling;
00135 
00136 BigReal*        ComputeNonbondedUtil::lambda_table = 0;
00137 
00138 Bool            ComputeNonbondedUtil::pairInteractionOn;
00139 Bool            ComputeNonbondedUtil::pairInteractionSelf;
00140 
00141 Bool            ComputeNonbondedUtil::pressureProfileOn;
00142 int             ComputeNonbondedUtil::pressureProfileSlabs;
00143 int             ComputeNonbondedUtil::pressureProfileAtomTypes;
00144 BigReal         ComputeNonbondedUtil::pressureProfileThickness;
00145 BigReal         ComputeNonbondedUtil::pressureProfileMin;
00146 
00147 Bool            ComputeNonbondedUtil::accelMDOn;
00148 
00149 Bool            ComputeNonbondedUtil::drudeNbthole;
00150 
00151 BigReal         ComputeNonbondedUtil::ewaldcof;
00152 BigReal         ComputeNonbondedUtil::pi_ewaldcof;
00153 
00154 int             ComputeNonbondedUtil::vdw_switch_mode;
00155 
00156 // Ported by JLai -- JE - Go
00157 Bool            ComputeNonbondedUtil::goGroPair;
00158 Bool            ComputeNonbondedUtil::goForcesOn;
00159 int             ComputeNonbondedUtil::goMethod; //6.3.11
00160 // End of port -- JLai
00161 
00162 void (*ComputeNonbondedUtil::calcPair)(nonbonded *);
00163 void (*ComputeNonbondedUtil::calcPairEnergy)(nonbonded *);
00164 void (*ComputeNonbondedUtil::calcSelf)(nonbonded *);
00165 void (*ComputeNonbondedUtil::calcSelfEnergy)(nonbonded *);
00166 
00167 void (*ComputeNonbondedUtil::calcFullPair)(nonbonded *);
00168 void (*ComputeNonbondedUtil::calcFullPairEnergy)(nonbonded *);
00169 void (*ComputeNonbondedUtil::calcFullSelf)(nonbonded *);
00170 void (*ComputeNonbondedUtil::calcFullSelfEnergy)(nonbonded *);
00171 
00172 void (*ComputeNonbondedUtil::calcMergePair)(nonbonded *);
00173 void (*ComputeNonbondedUtil::calcMergePairEnergy)(nonbonded *);
00174 void (*ComputeNonbondedUtil::calcMergeSelf)(nonbonded *);
00175 void (*ComputeNonbondedUtil::calcMergeSelfEnergy)(nonbonded *);
00176 
00177 void (*ComputeNonbondedUtil::calcSlowPair)(nonbonded *);
00178 void (*ComputeNonbondedUtil::calcSlowPairEnergy)(nonbonded *);
00179 void (*ComputeNonbondedUtil::calcSlowSelf)(nonbonded *);
00180 void (*ComputeNonbondedUtil::calcSlowSelfEnergy)(nonbonded *);
00181 
00182 // define splitting function
00183 #define SPLIT_NONE      1
00184 #define SPLIT_SHIFT     2
00185 #define SPLIT_C1        3
00186 #define SPLIT_XPLOR     4
00187 #define SPLIT_C2        5
00188 #define SPLIT_MARTINI   6
00189 
00190 void ComputeNonbondedUtil::submitReductionData(BigReal *data, SubmitReduction *reduction)
00191 {
00192   reduction->item(REDUCTION_EXCLUSION_CHECKSUM) += data[exclChecksumIndex];
00193   reduction->item(REDUCTION_PAIRLIST_WARNINGS) += data[pairlistWarningIndex];
00194   reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
00195   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += data[fullElectEnergyIndex];
00196   reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
00197   // Ported by JLai
00198   reduction->item(REDUCTION_GRO_LJ_ENERGY) += data[groLJEnergyIndex];
00199   reduction->item(REDUCTION_GRO_GAUSS_ENERGY) += data[groGaussEnergyIndex];
00200   reduction->item(REDUCTION_GO_NATIVE_ENERGY) += data[goNativeEnergyIndex];
00201   reduction->item(REDUCTION_GO_NONNATIVE_ENERGY) += data[goNonnativeEnergyIndex];
00202   // End of port -- JLai
00203 //fepb
00204   reduction->item(REDUCTION_ELECT_ENERGY_F) += data[electEnergyIndex_s];
00205   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += data[fullElectEnergyIndex_s];
00206   reduction->item(REDUCTION_LJ_ENERGY_F) += data[vdwEnergyIndex_s];
00207   reduction->item(REDUCTION_LJ_ENERGY_F_LEFT) += data[vdwEnergyIndex_s_Left];
00208 
00209   reduction->item(REDUCTION_ELECT_ENERGY_TI_1) += data[electEnergyIndex_ti_1];
00210   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1) += data[fullElectEnergyIndex_ti_1];
00211   reduction->item(REDUCTION_LJ_ENERGY_TI_1) += data[vdwEnergyIndex_ti_1];
00212   reduction->item(REDUCTION_ELECT_ENERGY_TI_2) += data[electEnergyIndex_ti_2];
00213   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2) += data[fullElectEnergyIndex_ti_2];
00214   reduction->item(REDUCTION_LJ_ENERGY_TI_2) += data[vdwEnergyIndex_ti_2];
00215 //fepe
00216   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00217   ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,fullElectVirialIndex);
00218   ADD_VECTOR(reduction,REDUCTION_PAIR_VDW_FORCE,data,pairVDWForceIndex);
00219   ADD_VECTOR(reduction,REDUCTION_PAIR_ELECT_FORCE,data,pairElectForceIndex);
00220   reduction->item(REDUCTION_COMPUTE_CHECKSUM) += 1.;
00221 }
00222 
00223 void ComputeNonbondedUtil::submitPressureProfileData(BigReal *data,
00224   SubmitReduction *reduction)
00225 {
00226   if (!reduction) return;
00227   int numAtomTypes = pressureProfileAtomTypes;
00228   // For ease of calculation we stored interactions between types
00229   // i and j in (ni+j).  For efficiency now we coalesce the
00230   // cross interactions so that just i<=j are stored.
00231   const int arraysize = 3*pressureProfileSlabs;
00232   size_t nelems = arraysize*(numAtomTypes*(numAtomTypes+1))/2;
00233   BigReal *arr = new BigReal[nelems];
00234   memset(arr, 0, nelems*sizeof(BigReal));
00235 
00236   int i, j;
00237   for (i=0; i<numAtomTypes; i++) {
00238     for (j=0; j<numAtomTypes; j++) {
00239       int ii=i;
00240       int jj=j;
00241       if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00242       const int reductionOffset = (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00243       for (int k=0; k<arraysize; k++) {
00244         arr[reductionOffset+k] += data[k];
00245       }
00246       data += arraysize;
00247     }
00248   }
00249   // copy into reduction
00250   reduction->add(nelems, arr);
00251   delete [] arr;
00252 }
00253   
00254 void ComputeNonbondedUtil::calc_error(nonbonded *) {
00255   NAMD_bug("Tried to call missing nonbonded compute routine.");
00256 }
00257   
00258 void ComputeNonbondedUtil::select(void)
00259 {
00260   if ( CkMyRank() ) return;
00261 
00262   // These defaults die cleanly if nothing appropriate is assigned.
00263   ComputeNonbondedUtil::calcPair = calc_error;
00264   ComputeNonbondedUtil::calcPairEnergy = calc_error;
00265   ComputeNonbondedUtil::calcSelf = calc_error;
00266   ComputeNonbondedUtil::calcSelfEnergy = calc_error;
00267   ComputeNonbondedUtil::calcFullPair = calc_error;
00268   ComputeNonbondedUtil::calcFullPairEnergy = calc_error;
00269   ComputeNonbondedUtil::calcFullSelf = calc_error;
00270   ComputeNonbondedUtil::calcFullSelfEnergy = calc_error;
00271   ComputeNonbondedUtil::calcMergePair = calc_error;
00272   ComputeNonbondedUtil::calcMergePairEnergy = calc_error;
00273   ComputeNonbondedUtil::calcMergeSelf = calc_error;
00274   ComputeNonbondedUtil::calcMergeSelfEnergy = calc_error;
00275   ComputeNonbondedUtil::calcSlowPair = calc_error;
00276   ComputeNonbondedUtil::calcSlowPairEnergy = calc_error;
00277   ComputeNonbondedUtil::calcSlowSelf = calc_error;
00278   ComputeNonbondedUtil::calcSlowSelfEnergy = calc_error;
00279 
00280   SimParameters * simParams = Node::Object()->simParameters;
00281   Parameters * params = Node::Object()->parameters;
00282 
00283   table_ener = params->table_ener;
00284   rowsize = params->rowsize;
00285   columnsize = params->columnsize;
00286 
00287   commOnly = simParams->commOnly;
00288   fixedAtomsOn = ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces );
00289 
00290   qmForcesOn = simParams->qmForcesOn ;
00291   
00292   cutoff = simParams->cutoff;
00293   cutoff2 = cutoff*cutoff;
00294   cutoff2_f = cutoff2;
00295 
00296 //fepb
00297   alchFepOn = simParams->alchFepOn;
00298   Fep_WCA_repuOn = simParams->alchFepWCARepuOn;
00299   Fep_WCA_dispOn = simParams->alchFepWCADispOn;
00300   Fep_ElecOn = simParams->alchFepElecOn;
00301   Fep_Wham = simParams->alchFepWhamOn;
00302   alchThermIntOn = simParams->alchThermIntOn;
00303   alchLambda2 = 0;
00304   lesOn = simParams->lesOn;
00305   lesScaling = lesFactor = 0;
00306   Bool tabulatedEnergies = simParams->tabulatedEnergies;
00307   alchVdwShiftCoeff = simParams->alchVdwShiftCoeff;
00308   vdwForceSwitching = simParams->vdwForceSwitching;
00309   LJcorrection = simParams->LJcorrection;
00310   WCA_rcut1 = simParams->alchFepWCArcut1;
00311   WCA_rcut2 = simParams->alchFepWCArcut2;
00312   WCA_rcut3 = simParams->alchFepWCArcut3;
00313 
00314   alchRepLambda = simParams->alchRepLambda;
00315   alchDispLambda = simParams->alchDispLambda;
00316   alchElecLambda = simParams->alchElecLambda;
00317 
00318   alchDecouple = simParams->alchDecouple;
00319 
00320   delete [] lambda_table;
00321   lambda_table = 0;
00322 
00323   pairInteractionOn = simParams->pairInteractionOn;
00324   pairInteractionSelf = simParams->pairInteractionSelf;
00325   pressureProfileOn = simParams->pressureProfileOn;
00326 
00327   // Ported by JLai -- Original JE - Go
00328   goGroPair = simParams->goGroPair;
00329   goForcesOn = simParams->goForcesOn;
00330   goMethod = simParams->goMethod; 
00331   // End of port
00332 
00333   accelMDOn = simParams->accelMDOn;
00334 
00335   drudeNbthole = simParams->drudeOn && (simParams->drudeNbtholeCut > 0.0);
00336 
00337   if ( drudeNbthole ) {
00338 #ifdef NAMD_CUDA
00339     NAMD_die("drudeNbthole is not supported in CUDA version");
00340 #endif
00341     if ( alchFepOn )
00342       NAMD_die("drudeNbthole is not supported with alchemical free-energy perturbation");
00343     if ( alchThermIntOn )
00344       NAMD_die("drudeNbthole is not supported with alchemical thermodynamic integration");
00345     if ( lesOn )
00346       NAMD_die("drudeNbthole is not supported with locally enhanced sampling");
00347     if ( pairInteractionOn )
00348       NAMD_die("drudeNbthole is not supported with pair interaction calculation");
00349     if ( pressureProfileOn )
00350       NAMD_die("drudeNbthole is not supported with pressure profile calculation");
00351   }
00352 
00353   if ( alchFepOn ) {
00354 #ifdef NAMD_CUDA
00355     NAMD_die("Alchemical free-energy perturbation is not supported in CUDA version");
00356 #endif
00357     alchLambda2 = simParams->alchLambda2;
00358     ComputeNonbondedUtil::calcPair = calc_pair_energy_fep;
00359     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_fep;
00360     ComputeNonbondedUtil::calcSelf = calc_self_energy_fep;
00361     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_fep;
00362     ComputeNonbondedUtil::calcFullPair = calc_pair_energy_fullelect_fep;
00363     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_fep;
00364     ComputeNonbondedUtil::calcFullSelf = calc_self_energy_fullelect_fep;
00365     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_fep;
00366     ComputeNonbondedUtil::calcMergePair = calc_pair_energy_merge_fullelect_fep;
00367     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_fep;
00368     ComputeNonbondedUtil::calcMergeSelf = calc_self_energy_merge_fullelect_fep;
00369     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_fep;
00370     ComputeNonbondedUtil::calcSlowPair = calc_pair_energy_slow_fullelect_fep;
00371     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_fep;
00372     ComputeNonbondedUtil::calcSlowSelf = calc_self_energy_slow_fullelect_fep;
00373     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_fep;
00374   }  else if ( alchThermIntOn ) {
00375 #ifdef NAMD_CUDA
00376     NAMD_die("Alchemical thermodynamic integration is not supported in CUDA version");
00377 #endif
00378     alchLambda2 = simParams->alchLambda2;
00379     ComputeNonbondedUtil::calcPair = calc_pair_ti;
00380     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_ti;
00381     ComputeNonbondedUtil::calcSelf = calc_self_ti;
00382     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_ti;
00383     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_ti;
00384     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_ti;
00385     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_ti;
00386     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_ti;
00387     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_ti;
00388     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_ti;
00389     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_ti;
00390     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_ti;
00391     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_ti;
00392     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_ti;
00393     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_ti;
00394     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_ti;
00395   } else if ( lesOn ) {
00396 #ifdef NAMD_CUDA
00397     NAMD_die("Locally enhanced sampling is not supported in CUDA version");
00398 #endif
00399     lesFactor = simParams->lesFactor;
00400     lesScaling = 1.0 / (double)lesFactor;
00401     lambda_table = new BigReal[(lesFactor+1)*(lesFactor+1)];
00402     for ( int ip=0; ip<=lesFactor; ++ip ) {
00403       for ( int jp=0; jp<=lesFactor; ++jp ) {
00404         BigReal lambda_pair = 1.0;
00405         if (ip || jp ) {
00406           if (ip && jp && ip != jp) {
00407             lambda_pair = 0.0;
00408           } else {
00409             lambda_pair = lesScaling;
00410           }
00411         }
00412         lambda_table[(lesFactor+1)*ip+jp] = lambda_pair;
00413       }
00414     }
00415     ComputeNonbondedUtil::calcPair = calc_pair_les;
00416     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_les;
00417     ComputeNonbondedUtil::calcSelf = calc_self_les;
00418     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_les;
00419     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_les;
00420     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_les;
00421     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_les;
00422     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_les;
00423     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_les;
00424     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_les;
00425     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_les;
00426     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_les;
00427     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_les;
00428     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_les;
00429     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_les;
00430     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_les;
00431   } else if ( pressureProfileOn) {
00432 #ifdef NAMD_CUDA
00433     NAMD_die("Pressure profile calculation is not supported in CUDA version");
00434 #endif
00435     pressureProfileSlabs = simParams->pressureProfileSlabs;
00436     pressureProfileAtomTypes = simParams->pressureProfileAtomTypes;
00437 
00438     ComputeNonbondedUtil::calcPair = calc_pair_pprof;
00439     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_pprof;
00440     ComputeNonbondedUtil::calcSelf = calc_self_pprof;
00441     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_pprof;
00442     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_pprof;
00443     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_pprof;
00444     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_pprof;
00445     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_pprof;
00446     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_pprof;
00447     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_pprof;
00448     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_pprof;
00449     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_pprof;
00450     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_pprof;
00451     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_pprof;
00452     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_pprof;
00453     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_pprof;
00454   } else if ( pairInteractionOn ) {
00455 #ifdef NAMD_CUDA
00456     NAMD_die("Pair interaction calculation is not supported in CUDA version");
00457 #endif
00458     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_int;
00459     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_int;
00460     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_int;
00461     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_int;
00462     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_int;
00463     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_int;
00464   } else if ( tabulatedEnergies ) {
00465 #ifdef NAMD_CUDA
00466     NAMD_die("Tabulated energies is not supported in CUDA version");
00467 #endif
00468     ComputeNonbondedUtil::calcPair = calc_pair_tabener;
00469     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_tabener;
00470     ComputeNonbondedUtil::calcSelf = calc_self_tabener;
00471     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_tabener;
00472     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_tabener;
00473     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_tabener;
00474     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_tabener;
00475     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_tabener;
00476     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_tabener;
00477     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_tabener;
00478     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_tabener;
00479     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_tabener;
00480     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_tabener;
00481     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_tabener;
00482     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_tabener;
00483     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_tabener;
00484   } else if ( goForcesOn ) {
00485 #ifdef NAMD_CUDA
00486     NAMD_die("Go forces is not supported in CUDA version");
00487 #endif
00488     ComputeNonbondedUtil::calcPair = calc_pair_go;
00489     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_go;
00490     ComputeNonbondedUtil::calcSelf = calc_self_go;
00491     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_go;
00492     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_go;
00493     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_go;
00494     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_go;
00495     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_go;
00496     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_go;
00497     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_go;
00498     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_go;
00499     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_go;
00500     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_go;
00501     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_go;
00502     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_go;
00503     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_go;
00504   } else {
00505     ComputeNonbondedUtil::calcPair = calc_pair;
00506     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy;
00507     ComputeNonbondedUtil::calcSelf = calc_self;
00508     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy;
00509     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect;
00510     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect;
00511     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect;
00512     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect;
00513     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect;
00514     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect;
00515     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect;
00516     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect;
00517     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect;
00518     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect;
00519     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect;
00520     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect;
00521   }
00522 
00523 //fepe
00524 
00525   dielectric_1 = 1.0/simParams->dielectric;
00526   if ( ! ljTable ) ljTable = new LJTable;
00527   mol = Node::Object()->molecule;
00528   scaling = simParams->nonbondedScaling;
00529   if ( simParams->exclude == SCALED14 )
00530   {
00531     scale14 = simParams->scale14;
00532   }
00533   else
00534   {
00535     scale14 = 1.;
00536   }
00537   if ( simParams->switchingActive )
00538   {
00539     switchOn = simParams->switchingDist;
00540     switchOn_1 = 1.0/switchOn;
00541     // d0 = 1.0/(cutoff-switchOn);
00542     switchOn2 = switchOn*switchOn;
00543     c0 = 1.0/(cutoff2-switchOn2);
00544 
00545     if ( simParams->vdwForceSwitching ) {
00546       double switchOn3 = switchOn * switchOn2;
00547       double cutoff3 = cutoff * cutoff2;
00548       double switchOn6 = switchOn3 * switchOn3;
00549       double cutoff6 = cutoff3 * cutoff3;
00550       v_vdwa_f = v_vdwa = -1. / ( switchOn6 * cutoff6 );
00551       v_vdwb_f = v_vdwb = -1. / ( switchOn3 * cutoff3 );
00552       k_vdwa_f = k_vdwa = cutoff6 / ( cutoff6 - switchOn6 );
00553       k_vdwb_f = k_vdwb = cutoff3 / ( cutoff3 - switchOn3 );
00554       cutoff_3_f = cutoff_3 = 1. / cutoff3;
00555       cutoff_6_f = cutoff_6 = 1. / cutoff6;
00556 
00557     } else if ( simParams->martiniSwitching ) { // switching fxn for Martini RBCG
00558 
00559       BigReal p6 = 6;
00560       BigReal A6 = p6 * ((p6+1)*switchOn-(p6+4)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,2));
00561       BigReal B6 = -p6 * ((p6+1)*switchOn-(p6+3)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,3));        
00562       BigReal C6 = 1.0/pow(cutoff,p6)-A6/3.0*pow(cutoff-switchOn,3)-B6/4.0*pow(cutoff-switchOn,4);
00563 
00564       BigReal p12 = 12;
00565       BigReal A12 = p12 * ((p12+1)*switchOn-(p12+4)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,2));
00566       BigReal B12 = -p12 * ((p12+1)*switchOn-(p12+3)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,3));
00567       BigReal C12 = 1.0/pow(cutoff,p12)-A12/3.0*pow(cutoff-switchOn,3)-B12/4.0*pow(cutoff-switchOn,4);
00568 
00569       A6_f =  A6;  B6_f  = B6;  C6_f =  C6;
00570       A12_f = A12; B12_f = B12; C12_f = C12;
00571       switchOn_f = switchOn;
00572 
00573     }
00574 
00575   }
00576   else
00577   {
00578     switchOn = cutoff;
00579     switchOn_1 = 1.0/switchOn;
00580     // d0 = 0.;  // avoid division by zero
00581     switchOn2 = switchOn*switchOn;
00582     c0 = 0.;  // avoid division by zero
00583   }
00584   c1 = c0*c0*c0;
00585   c3 = 3.0 * (cutoff2 - switchOn2);
00586   c5 = 0;
00587   c6 = 0;
00588   c7 = 0;
00589   c8 = 0;
00590 
00591   const int PMEOn = simParams->PMEOn;
00592   const int MSMOn = simParams->MSMOn;
00593   const int MSMSplit = simParams->MSMSplit;
00594 
00595   if ( PMEOn ) {
00596     ewaldcof = simParams->PMEEwaldCoefficient;
00597     BigReal TwoBySqrtPi = 1.12837916709551;
00598     pi_ewaldcof = TwoBySqrtPi * ewaldcof;
00599   }
00600 
00601   int splitType = SPLIT_NONE;
00602   if ( simParams->switchingActive ) splitType = SPLIT_SHIFT;
00603   if ( simParams->martiniSwitching ) splitType = SPLIT_MARTINI;
00604   if ( simParams->fullDirectOn || simParams->FMAOn || PMEOn || MSMOn ||
00605       simParams->FMMOn ) {
00606     switch ( simParams->longSplitting ) {
00607       case C2:
00608       splitType = SPLIT_C2;
00609       break;
00610 
00611       case C1:
00612       splitType = SPLIT_C1;
00613       break;
00614 
00615       case XPLOR:
00616       NAMD_die("Sorry, XPLOR splitting not supported.");
00617       break;
00618 
00619       case SHARP:
00620       NAMD_die("Sorry, SHARP splitting not supported.");
00621       break;
00622 
00623       default:
00624       NAMD_die("Unknown splitting type found!");
00625 
00626     }
00627   }
00628 
00629   BigReal r2_tol = 0.1;
00630   
00631   r2_delta = 1.0;
00632   r2_delta_exp = 0;
00633   while ( r2_delta > r2_tol ) { r2_delta /= 2.0; r2_delta_exp += 1; }
00634   r2_delta_1 = 1.0 / r2_delta;
00635 
00636   if ( ! CkMyPe() ) {
00637     iout << iINFO << "NONBONDED TABLE R-SQUARED SPACING: " <<
00638                                 r2_delta << "\n" << endi;
00639   }
00640 
00641   BigReal r2_tmp = 1.0;
00642   int cutoff2_exp = 0;
00643   while ( (cutoff2 + r2_delta) > r2_tmp ) { r2_tmp *= 2.0; cutoff2_exp += 1; }
00644 
00645   int i;
00646   int n = (r2_delta_exp + cutoff2_exp) * 64 + 1;
00647   table_length = n;
00648   #if defined(NAMD_MIC)
00649     int n_16 = (n + 15) & (~15);
00650   #endif
00651 
00652   if ( ! CkMyPe() ) {
00653     iout << iINFO << "NONBONDED TABLE SIZE: " <<
00654                                 n << " POINTS\n" << endi;
00655   }
00656 
00657   if ( table_alloc ) delete [] table_alloc;
00658   #if defined(NAMD_MIC)
00659     table_alloc = new BigReal[61*n_16+16];
00660     BigReal *table_align = table_alloc;
00661     while ( ((long)table_align) % 128 ) ++table_align;
00662     mic_table_base_ptr = table_align;
00663     mic_table_n = n;
00664     mic_table_n_16 = n_16;
00665     table_noshort = table_align;
00666     table_short = table_align + 16*n_16;
00667     slow_table = table_align + 32*n_16;
00668     fast_table = table_align + 36*n_16;
00669     scor_table = table_align + 40*n_16;
00670     corr_table = table_align + 44*n_16;
00671     full_table = table_align + 48*n_16;
00672     vdwa_table = table_align + 52*n_16;
00673     vdwb_table = table_align + 56*n_16;
00674     r2_table = table_align + 60*n_16;
00675   #else
00676   table_alloc = new BigReal[61*n+16];
00677   BigReal *table_align = table_alloc;
00678   while ( ((long)table_align) % 128 ) ++table_align;
00679   table_noshort = table_align;
00680   table_short = table_align + 16*n;
00681   slow_table = table_align + 32*n;
00682   fast_table = table_align + 36*n;
00683   scor_table = table_align + 40*n;
00684   corr_table = table_align + 44*n;
00685   full_table = table_align + 48*n;
00686   vdwa_table = table_align + 52*n;
00687   vdwb_table = table_align + 56*n;
00688   r2_table = table_align + 60*n;
00689   #endif
00690   BigReal *fast_i = fast_table + 4;
00691   BigReal *scor_i = scor_table + 4;
00692   BigReal *slow_i = slow_table + 4;
00693   BigReal *vdwa_i = vdwa_table + 4;
00694   BigReal *vdwb_i = vdwb_table + 4;
00695   BigReal *r2_i = r2_table;  *(r2_i++) = r2_delta;
00696   BigReal r2_limit = simParams->limitDist * simParams->limitDist;
00697   if ( r2_limit < r2_delta ) r2_limit = r2_delta;
00698   int r2_delta_i = 0;  // entry for r2 == r2_delta
00699 
00700 #ifdef NAMD_KNL
00701  if ( knl_table_alloc ) delete [] knl_table_alloc;
00702  knl_table_alloc = new float[10*KNL_TABLE_SIZE];
00703  knl_fast_ener_table = knl_table_alloc;
00704  knl_fast_grad_table = knl_table_alloc + KNL_TABLE_SIZE;
00705  knl_scor_ener_table = knl_table_alloc + 2*KNL_TABLE_SIZE;
00706  knl_scor_grad_table = knl_table_alloc + 3*KNL_TABLE_SIZE;
00707  knl_slow_ener_table = knl_table_alloc + 4*KNL_TABLE_SIZE;
00708  knl_slow_grad_table = knl_table_alloc + 5*KNL_TABLE_SIZE;
00709  knl_corr_ener_table = knl_table_alloc + 6*KNL_TABLE_SIZE;
00710  knl_corr_grad_table = knl_table_alloc + 7*KNL_TABLE_SIZE;
00711  knl_full_ener_table = knl_table_alloc + 8*KNL_TABLE_SIZE;
00712  knl_full_grad_table = knl_table_alloc + 9*KNL_TABLE_SIZE;
00713  knl_fast_ener_table[0] = 0.;
00714  knl_fast_grad_table[0] = 0.;
00715  knl_scor_ener_table[0] = 0.;
00716  knl_scor_grad_table[0] = 0.;
00717  knl_slow_ener_table[0] = 0.;
00718  knl_slow_grad_table[0] = 0.;
00719  knl_corr_ener_table[0] = 0.;
00720  knl_corr_grad_table[0] = 0.;
00721  knl_full_ener_table[0] = 0.;
00722  knl_full_grad_table[0] = 0.;
00723  for ( int knl_table = 0; knl_table < 2; ++knl_table ) {
00724   int nn = n;
00725   if ( knl_table ) {
00726     nn = KNL_TABLE_SIZE-1;
00727   }
00728   for ( i=1; i<nn; ++i ) {
00729 #else
00730   // fill in the table, fix up i==0 (r2==0) below
00731   for ( i=1; i<n; ++i ) {
00732 #endif
00733 
00734     const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00735     const BigReal r2_del = r2_base / 64.0;
00736     BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00737 
00738     BigReal r = sqrt(r2);
00739 
00740 #ifdef NAMD_KNL
00741     if ( knl_table ) {
00742       r = (double)(nn-1)/(double)(i);
00743       r2 = r*r;
00744     } else
00745 #endif
00746     if ( r2 <= r2_limit ) r2_delta_i = i;
00747 
00748     const BigReal r_1 = 1.0/r;
00749     const BigReal r_2 = 1.0/r2;
00750 
00751     // fast_ is defined as (full_ - slow_)
00752     // corr_ and fast_ are both zero at the cutoff, full_ is not
00753     // all three are approx 1/r at short distances
00754 
00755     // for actual interpolation, we use fast_ for fast forces and
00756     // scor_ = slow_ + corr_ - full_ and slow_ for slow forces
00757     // since these last two are of small magnitude
00758 
00759     BigReal fast_energy, fast_gradient;
00760     BigReal scor_energy, scor_gradient;
00761     BigReal slow_energy, slow_gradient;
00762 
00763     // corr_ is PME direct sum, or similar correction term
00764     // corr_energy is multiplied by r until later
00765     // corr_gradient is multiplied by -r^2 until later
00766     BigReal corr_energy, corr_gradient;
00767 
00768     
00769     if ( PMEOn ) {
00770       BigReal tmp_a = r * ewaldcof;
00771       BigReal tmp_b = erfc(tmp_a);
00772       corr_energy = tmp_b;
00773       corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
00774     } else if ( MSMOn ) {
00775       BigReal a_1 = 1.0/cutoff;
00776       BigReal r_a = r * a_1;
00777       BigReal g, dg;
00778       SPOLY(&g, &dg, r_a, MSMSplit);
00779       corr_energy = 1 - r_a * g;
00780       corr_gradient = 1 + r_a*r_a * dg;
00781     } else {
00782       corr_energy = corr_gradient = 0;
00783     }
00784 
00785     switch(splitType) {
00786       case SPLIT_NONE:
00787         fast_energy = 1.0/r;
00788         fast_gradient = -1.0/r2;
00789         scor_energy = scor_gradient = 0;
00790         slow_energy = slow_gradient = 0;
00791         break;
00792       case SPLIT_SHIFT: {
00793         BigReal shiftVal = r2/cutoff2 - 1.0;
00794         shiftVal *= shiftVal;
00795         BigReal dShiftVal = 2.0 * (r2/cutoff2 - 1.0) * 2.0*r/cutoff2;
00796         fast_energy = shiftVal/r;
00797         fast_gradient = dShiftVal/r - shiftVal/r2;
00798         scor_energy = scor_gradient = 0;
00799         slow_energy = slow_gradient = 0;
00800         } 
00801         break;
00802       case SPLIT_MARTINI: { 
00803         // in Martini, the Coulomb switching distance is zero
00804         const BigReal COUL_SWITCH = 0.;
00805         // Gromacs shifting function
00806         const BigReal p1 = 1.;
00807         BigReal A1 = p1 * ((p1+1)*COUL_SWITCH-(p1+4)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,2));
00808         BigReal B1 = -p1 * ((p1+1)*COUL_SWITCH-(p1+3)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,3));
00809         BigReal X1 = 1.0/pow(cutoff,p1)-A1/3.0*pow(cutoff-COUL_SWITCH,3)-B1/4.0*pow(cutoff-COUL_SWITCH,4);
00810         BigReal r12 = (r-COUL_SWITCH)*(r-COUL_SWITCH);
00811         BigReal r13 = (r-COUL_SWITCH)*(r-COUL_SWITCH)*(r-COUL_SWITCH);
00812         BigReal shiftVal = -(A1/3.0)*r13 - (B1/4.0)*r12*r12 - X1;
00813         BigReal dShiftVal = -A1*r12 - B1*r13;
00814         fast_energy = (1/r) + shiftVal;
00815         fast_gradient = -1/(r2) + dShiftVal;
00816         scor_energy = scor_gradient = 0;
00817         slow_energy = slow_gradient = 0;
00818         } 
00819         break;
00820       case SPLIT_C1:
00821         // calculate actual energy and gradient
00822         slow_energy = 0.5/cutoff * (3.0 - (r2/cutoff2));
00823         slow_gradient = -1.0/cutoff2 * (r/cutoff);
00824         // calculate scor from slow and corr
00825         scor_energy = slow_energy + (corr_energy - 1.0)/r;
00826         scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00827         // calculate fast from slow
00828         fast_energy = 1.0/r - slow_energy;
00829         fast_gradient = -1.0/r2 - slow_gradient;
00830         break;
00831       case SPLIT_C2:
00832         //
00833         // Quintic splitting function contributed by
00834         // Bruce Berne, Ruhong Zhou, and Joe Morrone
00835         //
00836         // calculate actual energy and gradient
00837         slow_energy = r2/(cutoff*cutoff2) * (6.0 * (r2/cutoff2)
00838             - 15.0*(r/cutoff) + 10.0);
00839         slow_gradient = r/(cutoff*cutoff2) * (24.0 * (r2/cutoff2)
00840             - 45.0 *(r/cutoff) + 20.0);
00841         // calculate scor from slow and corr
00842         scor_energy = slow_energy + (corr_energy - 1.0)/r;
00843         scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00844         // calculate fast from slow
00845         fast_energy = 1.0/r - slow_energy;
00846         fast_gradient = -1.0/r2 - slow_gradient;
00847         break;
00848     }
00849 
00850     // foo_gradient is calculated as ( d foo_energy / d r )
00851     // and now divided by 2r to get ( d foo_energy / d r2 )
00852 
00853     fast_gradient *= 0.5 * r_1;
00854     scor_gradient *= 0.5 * r_1;
00855     slow_gradient *= 0.5 * r_1;
00856 
00857     // let modf be 1 if excluded, 1-scale14 if modified, 0 otherwise,
00858     // add scor_ - modf * slow_ to slow terms and
00859     // add fast_ - modf * fast_ to fast terms.
00860 
00861     BigReal vdwa_energy, vdwa_gradient;
00862     BigReal vdwb_energy, vdwb_gradient;
00863 
00864     const BigReal r_6 = r_2*r_2*r_2;
00865     const BigReal r_12 = r_6*r_6;
00866 
00867     // Lennard-Jones switching function
00868   if ( simParams->vdwForceSwitching ) {  // switch force
00869     vdw_switch_mode = VDW_SWITCH_MODE_FORCE;
00870 
00871     // from Steinbach & Brooks, JCC 15, pgs 667-683, 1994, eqns 10-13
00872     if ( r2 > switchOn2 ) {
00873       BigReal tmpa = r_6 - cutoff_6;
00874       vdwa_energy = k_vdwa * tmpa * tmpa - (LJcorrection ? v_vdwa : 0.0);
00875       BigReal tmpb = r_1 * r_2 - cutoff_3;
00876       vdwb_energy = k_vdwb * tmpb * tmpb - (LJcorrection ? v_vdwb : 0.0);
00877       vdwa_gradient = -6.0 * k_vdwa * tmpa * r_2 * r_6;
00878       vdwb_gradient = -3.0 * k_vdwb * tmpb * r_2 * r_2 * r_1;
00879     } else {
00880       vdwa_energy = r_12 + (LJcorrection ? 0.0 : v_vdwa);
00881       vdwb_energy = r_6 + (LJcorrection ? 0.0 : v_vdwb);
00882       vdwa_gradient = -6.0 * r_2 * r_12;
00883       vdwb_gradient = -3.0 * r_2 * r_6;
00884     }
00885   } else if ( simParams->martiniSwitching ) { // switching fxn for Martini RBCG
00886     vdw_switch_mode = VDW_SWITCH_MODE_MARTINI;
00887 
00888     BigReal r12 = (r-switchOn)*(r-switchOn);        BigReal r13 = (r-switchOn)*(r-switchOn)*(r-switchOn);
00889 
00890     BigReal p6 = 6;
00891     BigReal A6 = p6 * ((p6+1)*switchOn-(p6+4)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,2));
00892     BigReal B6 = -p6 * ((p6+1)*switchOn-(p6+3)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,3));        
00893     BigReal C6 = 1.0/pow(cutoff,p6)-A6/3.0*pow(cutoff-switchOn,3)-B6/4.0*pow(cutoff-switchOn,4);
00894 
00895     BigReal p12 = 12;
00896     BigReal A12 = p12 * ((p12+1)*switchOn-(p12+4)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,2));
00897     BigReal B12 = -p12 * ((p12+1)*switchOn-(p12+3)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,3));
00898     BigReal C12 = 1.0/pow(cutoff,p12)-A12/3.0*pow(cutoff-switchOn,3)-B12/4.0*pow(cutoff-switchOn,4);
00899 
00900     BigReal LJshifttempA = -(A12/3)*r13 - (B12/4)*r12*r12 - C12;
00901     BigReal LJshifttempB = -(A6/3)*r13 - (B6/4)*r12*r12 - C6;
00902     const BigReal shiftValA =         // used for Lennard-Jones
00903                         ( r2 > switchOn2 ? LJshifttempA : -C12);
00904     const BigReal shiftValB =         // used for Lennard-Jones
00905                         ( r2 > switchOn2 ? LJshifttempB : -C6);
00906 
00907     BigReal LJdshifttempA = -A12*r12 - B12*r13;
00908     BigReal LJdshifttempB = -A6*r12 - B6*r13;
00909     const BigReal dshiftValA =         // used for Lennard-Jones
00910                         ( r2 > switchOn2 ? LJdshifttempA*0.5*r_1 : 0 );
00911     const BigReal dshiftValB =         // used for Lennard-Jones
00912                         ( r2 > switchOn2 ? LJdshifttempB*0.5*r_1 : 0 );
00913 
00914 
00915 
00916 
00917     //have not addressed r > cutoff
00918 
00919     //  dshiftValA*= 0.5*r_1;
00920     //  dshiftValB*= 0.5*r_1;
00921 
00922     vdwa_energy = r_12 + shiftValA;
00923     vdwb_energy = r_6 + shiftValB;
00924    
00925     vdwa_gradient = -6/pow(r,14) + dshiftValA ;
00926     vdwb_gradient = -3/pow(r,8) + dshiftValB;
00927 
00928   } else {  // switch energy
00929     vdw_switch_mode = VDW_SWITCH_MODE_ENERGY;
00930 
00931     const BigReal c2 = cutoff2-r2;
00932     const BigReal c4 = c2*(c3-2.0*c2);
00933     const BigReal switchVal =         // used for Lennard-Jones
00934                         ( r2 > switchOn2 ? c2*c4*c1 : 1.0 );
00935     const BigReal dSwitchVal =        // d switchVal / d r2
00936                         ( r2 > switchOn2 ? 2*c1*(c2*c2-c4) : 0.0 );
00937 
00938     vdwa_energy = switchVal * r_12;
00939     vdwb_energy = switchVal * r_6;
00940 
00941     vdwa_gradient = ( dSwitchVal - 6.0 * switchVal * r_2 ) * r_12;
00942     vdwb_gradient = ( dSwitchVal - 3.0 * switchVal * r_2 ) * r_6;
00943   }
00944 
00945 
00946 #ifdef NAMD_KNL
00947    if ( knl_table ) {
00948     knl_fast_ener_table[i] = -1.*fast_energy;
00949     knl_fast_grad_table[i] = -2.*fast_gradient;
00950     knl_scor_ener_table[i] = -1.*scor_energy;
00951     knl_scor_grad_table[i] = -2.*scor_gradient;
00952     knl_slow_ener_table[i] = -1.*slow_energy;
00953     knl_slow_grad_table[i] = -2.*slow_gradient;
00954     knl_corr_ener_table[i] = -1.*(fast_energy + scor_energy);
00955     knl_corr_grad_table[i] = -2.*(fast_gradient + scor_gradient);
00956     knl_full_ener_table[i] = -1.*(fast_energy + slow_energy);
00957     knl_full_grad_table[i] = -2.*(fast_gradient + slow_gradient);
00958     if ( i == nn-1 ) {
00959       knl_fast_ener_table[nn] = knl_fast_ener_table[i];
00960       knl_fast_grad_table[nn] = knl_fast_grad_table[i];
00961       knl_scor_ener_table[nn] = knl_scor_ener_table[i];
00962       knl_scor_grad_table[nn] = knl_scor_grad_table[i];
00963       knl_slow_ener_table[nn] = knl_slow_ener_table[i];
00964       knl_slow_grad_table[nn] = knl_slow_grad_table[i];
00965       knl_corr_ener_table[nn] = knl_corr_ener_table[i];
00966       knl_corr_grad_table[nn] = knl_corr_grad_table[i];
00967       knl_full_ener_table[nn] = knl_full_ener_table[i];
00968       knl_full_grad_table[nn] = knl_full_grad_table[i];
00969     }
00970    } else {
00971 #endif
00972     *(fast_i++) = fast_energy;
00973     *(fast_i++) = fast_gradient;
00974     *(fast_i++) = 0;
00975     *(fast_i++) = 0;
00976     *(scor_i++) = scor_energy;
00977     *(scor_i++) = scor_gradient;
00978     *(scor_i++) = 0;
00979     *(scor_i++) = 0;
00980     *(slow_i++) = slow_energy;
00981     *(slow_i++) = slow_gradient;
00982     *(slow_i++) = 0;
00983     *(slow_i++) = 0;
00984     *(vdwa_i++) = vdwa_energy;
00985     *(vdwa_i++) = vdwa_gradient;
00986     *(vdwa_i++) = 0;
00987     *(vdwa_i++) = 0;
00988     *(vdwb_i++) = vdwb_energy;
00989     *(vdwb_i++) = vdwb_gradient;
00990     *(vdwb_i++) = 0;
00991     *(vdwb_i++) = 0;
00992     *(r2_i++) = r2 + r2_delta;
00993 #ifdef NAMD_KNL
00994    }
00995 #endif
00996 
00997   }
00998 #ifdef NAMD_KNL
00999  } // knl_table loop
01000 #endif
01001 
01002   if ( ! r2_delta_i ) {
01003     NAMD_bug("Failed to find table entry for r2 == r2_limit\n");
01004   }
01005   if ( r2_table[r2_delta_i] > r2_limit + r2_delta ) {
01006     NAMD_bug("Found bad table entry for r2 == r2_limit\n");
01007   }
01008 
01009   int j;
01010   const char *table_name = "XXXX";
01011   int smooth_short = 0;
01012   for ( j=0; j<5; ++j ) {
01013     BigReal *t0 = 0;
01014     switch (j) {
01015       case 0: 
01016         t0 = fast_table;
01017         table_name = "FAST";
01018         smooth_short = 1;
01019       break;
01020       case 1: 
01021         t0 = scor_table;
01022         table_name = "SCOR";
01023         smooth_short = 0;
01024       break;
01025       case 2: 
01026         t0 = slow_table;
01027         table_name = "SLOW";
01028         smooth_short = 0;
01029       break;
01030       case 3: 
01031         t0 = vdwa_table;
01032         table_name = "VDWA";
01033         smooth_short = 1;
01034       break;
01035       case 4: 
01036         t0 = vdwb_table;
01037         table_name = "VDWB";
01038         smooth_short = 1;
01039       break;
01040     }
01041     // patch up data for i=0
01042     t0[0] = t0[4] - t0[5] * ( r2_delta / 64.0 );  // energy
01043     t0[1] = t0[5];  // gradient
01044     t0[2] = 0;
01045     t0[3] = 0;
01046     if ( smooth_short ) {
01047       BigReal energy0 = t0[4*r2_delta_i];
01048       BigReal gradient0 = t0[4*r2_delta_i+1];
01049       BigReal r20 = r2_table[r2_delta_i];
01050       t0[0] = energy0 - gradient0 * (r20 - r2_table[0]);  // energy
01051       t0[1] = gradient0;  // gradient
01052     }
01053     BigReal *t;
01054     for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
01055       BigReal x = ( r2_delta * ( 1 << (i/64) ) ) / 64.0;
01056       if ( r2_table[i+1] != r2_table[i] + x ) {
01057         NAMD_bug("Bad table delta calculation.\n");
01058       }
01059       if ( smooth_short && i+1 < r2_delta_i ) {
01060         BigReal energy0 = t0[4*r2_delta_i];
01061         BigReal gradient0 = t0[4*r2_delta_i+1];
01062         BigReal r20 = r2_table[r2_delta_i];
01063         t[4] = energy0 - gradient0 * (r20 - r2_table[i+1]);  // energy
01064         t[5] = gradient0;  // gradient
01065       }
01066       BigReal v1 = t[0];
01067       BigReal g1 = t[1];
01068       BigReal v2 = t[4];
01069       BigReal g2 = t[5];
01070       // explicit formulas for v1 + g1 x + c x^2 + d x^3
01071       BigReal c = ( 3.0 * (v2 - v1) - x * (2.0 * g1 + g2) ) / ( x * x );
01072       BigReal d = ( -2.0 * (v2 - v1) + x * (g1 + g2) ) / ( x * x * x );
01073       // since v2 - v1 is imprecise, we refine c and d numerically
01074       // important because we need accurate forces (more than energies!)
01075       for ( int k=0; k < 2; ++k ) {
01076         BigReal dv = (v1 - v2) + ( ( d * x + c ) * x + g1 ) * x;
01077         BigReal dg = (g1 - g2) + ( 3.0 * d * x + 2.0 * c ) * x;
01078         c -= ( 3.0 * dv - x * dg ) / ( x * x );
01079         d -= ( -2.0 * dv + x * dg ) / ( x * x * x );
01080       }
01081       // store in the array;
01082       t[2] = c;  t[3] = d;
01083     }
01084 
01085     if ( ! CkMyPe() ) {
01086     BigReal dvmax = 0;
01087     BigReal dgmax = 0;
01088     BigReal dvmax_r = 0;
01089     BigReal dgmax_r = 0;
01090     BigReal fdvmax = 0;
01091     BigReal fdgmax = 0;
01092     BigReal fdvmax_r = 0;
01093     BigReal fdgmax_r = 0;
01094     BigReal dgcdamax = 0;
01095     BigReal dgcdimax = 0;
01096     BigReal dgcaimax = 0;
01097     BigReal dgcdamax_r = 0;
01098     BigReal dgcdimax_r = 0;
01099     BigReal dgcaimax_r = 0;
01100     BigReal fdgcdamax = 0;
01101     BigReal fdgcdimax = 0;
01102     BigReal fdgcaimax = 0;
01103     BigReal fdgcdamax_r = 0;
01104     BigReal fdgcdimax_r = 0;
01105     BigReal fdgcaimax_r = 0;
01106     BigReal gcm = fabs(t0[1]);  // gradient magnitude running average
01107     for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
01108       const BigReal r2_base = r2_delta * ( 1 << (i/64) );
01109       const BigReal r2_del = r2_base / 64.0;
01110       const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
01111       const BigReal r = sqrt(r2);
01112       if ( r > cutoff ) break;
01113       BigReal x = r2_del;
01114       BigReal dv = ( ( t[3] * x + t[2] ) * x + t[1] ) * x + t[0] - t[4];
01115       BigReal dg = ( 3.0 * t[3] * x + 2.0 * t[2] ) * x + t[1] - t[5];
01116       if ( t[4] != 0. && fabs(dv/t[4]) > fdvmax ) {
01117         fdvmax = fabs(dv/t[4]); fdvmax_r = r;
01118       }
01119       if ( fabs(dv) > dvmax ) {
01120         dvmax = fabs(dv); dvmax_r = r;
01121       }
01122       if ( t[5] != 0. && fabs(dg/t[5]) > fdgmax ) {
01123         fdgmax = fabs(dg/t[5]); fdgmax_r = r;
01124       }
01125       if ( fabs(dg) > dgmax ) {
01126         dgmax = fabs(dg); dgmax_r = r;
01127       }
01128       BigReal gcd = (t[4] - t[0]) / x;  // centered difference gradient
01129       BigReal gcd_prec = (fabs(t[0]) + fabs(t[4])) * 1.e-15 / x;  // roundoff
01130       gcm = 0.9 * gcm + 0.1 * fabs(t[5]);  // magnitude running average
01131       BigReal gca = 0.5  * (t[1] + t[5]);  // centered average gradient
01132       BigReal gci = ( 0.75 * t[3] * x + t[2] ) * x + t[1];  // interpolated
01133       BigReal rc = sqrt(r2 + 0.5 * x);
01134       BigReal dgcda = gcd - gca;
01135       if ( dgcda != 0. && fabs(dgcda) < gcd_prec ) {
01136         // CkPrintf("ERROR %g < PREC %g AT %g AVG VAL %g\n", dgcda, gcd_prec, rc, gca);
01137         dgcda = 0.;
01138       }
01139       BigReal dgcdi = gcd - gci;
01140       if ( dgcdi != 0. && fabs(dgcdi) < gcd_prec ) {
01141         // CkPrintf("ERROR %g < PREC %g AT %g INT VAL %g\n", dgcdi, gcd_prec, rc, gci);
01142         dgcdi = 0.;
01143       }
01144       BigReal dgcai = gca - gci;
01145       if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcda/gcm) > fdgcdamax ) {
01146         fdgcdamax = fabs(dgcda/gcm); fdgcdamax_r = rc;
01147       }
01148       if ( fabs(dgcda) > fdgcdamax ) {
01149         dgcdamax = fabs(dgcda); dgcdamax_r = rc;
01150       }
01151       if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcdi/gcm) > fdgcdimax ) {
01152         fdgcdimax = fabs(dgcdi/gcm); fdgcdimax_r = rc;
01153       }
01154       if ( fabs(dgcdi) > fdgcdimax ) {
01155         dgcdimax = fabs(dgcdi); dgcdimax_r = rc;
01156       }
01157       if ( t[1]*t[5] > 0. && gcm != 0. && fabs(dgcai/gcm) > fdgcaimax ) {
01158         fdgcaimax = fabs(dgcai/gcm); fdgcaimax_r = rc;
01159       }
01160       if ( fabs(dgcai) > fdgcaimax ) {
01161         dgcaimax = fabs(dgcai); dgcaimax_r = rc;
01162       }
01163 #if 0
01164       CkPrintf("TABLE %s %g %g %g %g\n",table_name,rc,dgcda/gcm,dgcda,gci);
01165       if (dv != 0.) CkPrintf("TABLE %d ENERGY ERROR %g AT %g (%d)\n",j,dv,r,i);
01166       if (dg != 0.) CkPrintf("TABLE %d FORCE ERROR %g AT %g (%d)\n",j,dg,r,i);
01167 #endif
01168     }
01169     if ( dvmax != 0.0 ) {
01170       iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
01171         " TABLE ENERGY: " << dvmax << " AT " << dvmax_r << "\n" << endi;
01172     }
01173     if ( fdvmax != 0.0 ) {
01174       iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
01175         " TABLE ENERGY: " << fdvmax << " AT " << fdvmax_r << "\n" << endi;
01176     }
01177     if ( dgmax != 0.0 ) {
01178       iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
01179         " TABLE FORCE: " << dgmax << " AT " << dgmax_r << "\n" << endi;
01180     }
01181     if ( fdgmax != 0.0 ) {
01182       iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
01183         " TABLE FORCE: " << fdgmax << " AT " << fdgmax_r << "\n" << endi;
01184     }
01185     if (fdgcdamax != 0.0 ) {
01186       iout << iINFO << "INCONSISTENCY IN " << table_name <<
01187         " TABLE ENERGY VS FORCE: " << fdgcdamax << " AT " << fdgcdamax_r << "\n" << endi;
01188       if ( fdgcdamax > 0.1 ) {
01189         iout << iERROR << "\n";
01190         iout << iERROR << "CALCULATED " << table_name <<
01191           " FORCE MAY NOT MATCH ENERGY! POSSIBLE BUG!\n";
01192         iout << iERROR << "\n";
01193       }
01194     }
01195     if (0 && fdgcdimax != 0.0 ) {
01196       iout << iINFO << "INCONSISTENCY IN " << table_name <<
01197         " TABLE ENERGY VS FORCE: " << fdgcdimax << " AT " << fdgcdimax_r << "\n" << endi;
01198     }
01199     if ( 0 && fdgcaimax != 0.0 ) {
01200       iout << iINFO << "INCONSISTENCY IN " << table_name <<
01201         " TABLE AVG VS INT FORCE: " << fdgcaimax << " AT " << fdgcaimax_r << "\n" << endi;
01202     }
01203     }
01204 
01205   }
01206 
01207   for ( i=0; i<4*n; ++i ) {
01208     corr_table[i] = fast_table[i] + scor_table[i];
01209     full_table[i] = fast_table[i] + slow_table[i];
01210   }
01211 
01212 #if 0  
01213   for ( i=0; i<n; ++i ) {
01214    for ( int j=0; j<4; ++j ) {
01215     table_short[16*i+6-2*j] = table_noshort[16*i+6-2*j] = vdwa_table[4*i+j];
01216     table_short[16*i+7-2*j] = table_noshort[16*i+7-2*j] = vdwb_table[4*i+j];
01217     table_short[16*i+8+3-j] = fast_table[4*i+j];
01218     table_short[16*i+12+3-j] = scor_table[4*i+j];
01219     table_noshort[16*i+8+3-j] = corr_table[4*i+j];
01220     table_noshort[16*i+12+3-j] = full_table[4*i+j];
01221    }
01222   }
01223 #endif 
01224 
01225   for ( i=0; i<n; ++i ) {
01226     table_short[16*i+0] = table_noshort[16*i+0] = -6.*vdwa_table[4*i+3];
01227     table_short[16*i+1] = table_noshort[16*i+1] = -4.*vdwa_table[4*i+2];
01228     table_short[16*i+2] = table_noshort[16*i+2] = -2.*vdwa_table[4*i+1];
01229     table_short[16*i+3] = table_noshort[16*i+3] = -1.*vdwa_table[4*i+0];
01230     
01231     table_short[16*i+4] = table_noshort[16*i+4] = -6.*vdwb_table[4*i+3];
01232     table_short[16*i+5] = table_noshort[16*i+5] = -4.*vdwb_table[4*i+2];
01233     table_short[16*i+6] = table_noshort[16*i+6] = -2.*vdwb_table[4*i+1];
01234     table_short[16*i+7] = table_noshort[16*i+7] = -1.*vdwb_table[4*i+0];
01235     
01236     table_short[16*i+8]  = -6.*fast_table[4*i+3];
01237     table_short[16*i+9]  = -4.*fast_table[4*i+2];
01238     table_short[16*i+10] = -2.*fast_table[4*i+1];
01239     table_short[16*i+11] = -1.*fast_table[4*i+0];
01240 
01241     table_noshort[16*i+8]  = -6.*corr_table[4*i+3];
01242     table_noshort[16*i+9]  = -4.*corr_table[4*i+2];
01243     table_noshort[16*i+10] = -2.*corr_table[4*i+1];
01244     table_noshort[16*i+11] = -1.*corr_table[4*i+0];
01245 
01246     table_short[16*i+12] = -6.*scor_table[4*i+3];
01247     table_short[16*i+13] = -4.*scor_table[4*i+2];
01248     table_short[16*i+14] = -2.*scor_table[4*i+1];
01249     table_short[16*i+15] = -1.*scor_table[4*i+0];
01250 
01251     table_noshort[16*i+12] = -6.*full_table[4*i+3];
01252     table_noshort[16*i+13] = -4.*full_table[4*i+2];
01253     table_noshort[16*i+14] = -2.*full_table[4*i+1];
01254     table_noshort[16*i+15] = -1.*full_table[4*i+0];
01255   }
01256 
01257 #if 0
01258   char fname[100];
01259   sprintf(fname,"/tmp/namd.table.pe%d.dat",CkMyPe());
01260   FILE *f = fopen(fname,"w");
01261   for ( i=0; i<(n-1); ++i ) {
01262     const BigReal r2_base = r2_delta * ( 1 << (i/64) );
01263     const BigReal r2_del = r2_base / 64.0;
01264     const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
01265     BigReal *t;
01266     if ( r2 + r2_delta != r2_table[i] ) fprintf(f,"r2 error! ");
01267     fprintf(f,"%g",r2);
01268     t = fast_table + 4*i;
01269     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01270     t = scor_table + 4*i;
01271     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01272     t = slow_table + 4*i;
01273     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01274     t = corr_table + 4*i;
01275     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01276     t = full_table + 4*i;
01277     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01278     t = vdwa_table + 4*i;
01279     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01280     t = vdwb_table + 4*i;
01281     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
01282     fprintf(f,"\n");
01283   }
01284   fclose(f);
01285 #endif
01286 
01287   //Flip slow table to match table_four_i
01288   for ( i=0; i<n; ++i ) {
01289     BigReal tmp0, tmp1, tmp2, tmp3;
01290     tmp0 = slow_table [i*4 + 0];
01291     tmp1 = slow_table [i*4 + 1];
01292     tmp2 = slow_table [i*4 + 2];
01293     tmp3 = slow_table [i*4 + 3];
01294 
01295     slow_table [i*4 + 0] = tmp3;
01296     slow_table [i*4 + 1] = tmp2;
01297     slow_table [i*4 + 2] = tmp1;
01298     slow_table [i*4 + 3] = tmp0;
01299   }
01300 
01301 #ifdef NAMD_CUDA
01302   if (!simParams->useCUDA2) {
01303     send_build_cuda_force_table();
01304   }
01305 #endif
01306 
01307   #ifdef NAMD_MIC
01308     send_build_mic_force_table();
01309   #endif
01310 }
01311 

Generated on Mon Nov 20 01:17:11 2017 for NAMD by  doxygen 1.4.7