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

Generated on Mon Jul 16 01:17:12 2018 for NAMD by  doxygen 1.4.7