Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

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 <stdio.h>
00026 
00027 #ifdef NAMD_CUDA
00028   void build_cuda_force_table();
00029 #endif
00030 
00031 Bool            ComputeNonbondedUtil::commOnly;
00032 Bool            ComputeNonbondedUtil::fixedAtomsOn;
00033 BigReal         ComputeNonbondedUtil::cutoff;
00034 BigReal         ComputeNonbondedUtil::cutoff2;
00035 BigReal         ComputeNonbondedUtil::dielectric_1;
00036 const LJTable*  ComputeNonbondedUtil::ljTable = 0;
00037 const Molecule* ComputeNonbondedUtil::mol;
00038 BigReal         ComputeNonbondedUtil::r2_delta;
00039 BigReal         ComputeNonbondedUtil::r2_delta_1;
00040 int             ComputeNonbondedUtil::r2_delta_exp;
00041 int             ComputeNonbondedUtil::rowsize;
00042 int             ComputeNonbondedUtil::columnsize;
00043 BigReal*        ComputeNonbondedUtil::table_alloc = 0;
00044 BigReal*                ComputeNonbondedUtil::table_ener = 0;
00045 BigReal*        ComputeNonbondedUtil::table_short;
00046 BigReal*        ComputeNonbondedUtil::table_noshort;
00047 BigReal*        ComputeNonbondedUtil::fast_table;
00048 BigReal*        ComputeNonbondedUtil::scor_table;
00049 BigReal*        ComputeNonbondedUtil::slow_table;
00050 BigReal*        ComputeNonbondedUtil::corr_table;
00051 BigReal*        ComputeNonbondedUtil::full_table;
00052 BigReal*        ComputeNonbondedUtil::vdwa_table;
00053 BigReal*        ComputeNonbondedUtil::vdwb_table;
00054 BigReal*        ComputeNonbondedUtil::r2_table;
00055 BigReal         ComputeNonbondedUtil::scaling;
00056 BigReal         ComputeNonbondedUtil::scale14;
00057 BigReal         ComputeNonbondedUtil::switchOn;
00058 BigReal         ComputeNonbondedUtil::switchOn_1;
00059 BigReal         ComputeNonbondedUtil::switchOn2;
00060 BigReal         ComputeNonbondedUtil::c0;
00061 BigReal         ComputeNonbondedUtil::c1;
00062 BigReal         ComputeNonbondedUtil::c3;
00063 BigReal         ComputeNonbondedUtil::c5;
00064 BigReal         ComputeNonbondedUtil::c6;
00065 BigReal         ComputeNonbondedUtil::c7;
00066 BigReal         ComputeNonbondedUtil::c8;
00067 // BigReal         ComputeNonbondedUtil::d0;
00068 // fepb
00069 Bool      ComputeNonbondedUtil::alchFepOn;
00070 Bool      ComputeNonbondedUtil::alchThermIntOn;
00071 BigReal   ComputeNonbondedUtil::alchLambda;
00072 BigReal   ComputeNonbondedUtil::alchLambda2;
00073 BigReal   ComputeNonbondedUtil::alchVdwShiftCoeff;
00074 BigReal   ComputeNonbondedUtil::alchVdwLambdaEnd;
00075 BigReal   ComputeNonbondedUtil::alchElecLambdaStart;
00076 Bool      ComputeNonbondedUtil::alchDecouple;
00077 //fepe
00078 Bool      ComputeNonbondedUtil::lesOn;
00079 int       ComputeNonbondedUtil::lesFactor;
00080 BigReal   ComputeNonbondedUtil::lesScaling;
00081 
00082 BigReal*        ComputeNonbondedUtil::lambda_table = 0;
00083 
00084 Bool            ComputeNonbondedUtil::pairInteractionOn;
00085 Bool            ComputeNonbondedUtil::pairInteractionSelf;
00086 
00087 Bool            ComputeNonbondedUtil::pressureProfileOn;
00088 int             ComputeNonbondedUtil::pressureProfileSlabs;
00089 int             ComputeNonbondedUtil::pressureProfileAtomTypes;
00090 BigReal         ComputeNonbondedUtil::pressureProfileThickness;
00091 BigReal         ComputeNonbondedUtil::pressureProfileMin;
00092 
00093 BigReal         ComputeNonbondedUtil::ewaldcof;
00094 BigReal         ComputeNonbondedUtil::pi_ewaldcof;
00095 
00096 void (*ComputeNonbondedUtil::calcPair)(nonbonded *);
00097 void (*ComputeNonbondedUtil::calcPairEnergy)(nonbonded *);
00098 void (*ComputeNonbondedUtil::calcSelf)(nonbonded *);
00099 void (*ComputeNonbondedUtil::calcSelfEnergy)(nonbonded *);
00100 
00101 void (*ComputeNonbondedUtil::calcFullPair)(nonbonded *);
00102 void (*ComputeNonbondedUtil::calcFullPairEnergy)(nonbonded *);
00103 void (*ComputeNonbondedUtil::calcFullSelf)(nonbonded *);
00104 void (*ComputeNonbondedUtil::calcFullSelfEnergy)(nonbonded *);
00105 
00106 void (*ComputeNonbondedUtil::calcMergePair)(nonbonded *);
00107 void (*ComputeNonbondedUtil::calcMergePairEnergy)(nonbonded *);
00108 void (*ComputeNonbondedUtil::calcMergeSelf)(nonbonded *);
00109 void (*ComputeNonbondedUtil::calcMergeSelfEnergy)(nonbonded *);
00110 
00111 void (*ComputeNonbondedUtil::calcSlowPair)(nonbonded *);
00112 void (*ComputeNonbondedUtil::calcSlowPairEnergy)(nonbonded *);
00113 void (*ComputeNonbondedUtil::calcSlowSelf)(nonbonded *);
00114 void (*ComputeNonbondedUtil::calcSlowSelfEnergy)(nonbonded *);
00115 
00116 // define splitting function
00117 #define SPLIT_NONE      1
00118 #define SPLIT_SHIFT     2
00119 #define SPLIT_C1        3
00120 #define SPLIT_XPLOR     4
00121 #define SPLIT_C2        5
00122 
00123 void ComputeNonbondedUtil::submitReductionData(BigReal *data, SubmitReduction *reduction)
00124 {
00125   reduction->item(REDUCTION_EXCLUSION_CHECKSUM) += data[exclChecksumIndex];
00126   reduction->item(REDUCTION_PAIRLIST_WARNINGS) += data[pairlistWarningIndex];
00127   reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
00128   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += data[fullElectEnergyIndex];
00129   reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
00130 //fepb
00131   reduction->item(REDUCTION_ELECT_ENERGY_F) += data[electEnergyIndex_s];
00132   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += data[fullElectEnergyIndex_s];
00133   reduction->item(REDUCTION_LJ_ENERGY_F) += data[vdwEnergyIndex_s];
00134 
00135   reduction->item(REDUCTION_ELECT_ENERGY_TI_1) += data[electEnergyIndex_ti_1];
00136   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1) += data[fullElectEnergyIndex_ti_1];
00137   reduction->item(REDUCTION_LJ_ENERGY_TI_1) += data[vdwEnergyIndex_ti_1];
00138   reduction->item(REDUCTION_ELECT_ENERGY_TI_2) += data[electEnergyIndex_ti_2];
00139   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2) += data[fullElectEnergyIndex_ti_2];
00140   reduction->item(REDUCTION_LJ_ENERGY_TI_2) += data[vdwEnergyIndex_ti_2];
00141 //fepe
00142   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00143   ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,fullElectVirialIndex);
00144   ADD_VECTOR(reduction,REDUCTION_PAIR_VDW_FORCE,data,pairVDWForceIndex);
00145   ADD_VECTOR(reduction,REDUCTION_PAIR_ELECT_FORCE,data,pairElectForceIndex);
00146   reduction->item(REDUCTION_COMPUTE_CHECKSUM) += 1.;
00147 }
00148 
00149 void ComputeNonbondedUtil::submitPressureProfileData(BigReal *data,
00150   SubmitReduction *reduction)
00151 {
00152   if (!reduction) return;
00153   int numAtomTypes = pressureProfileAtomTypes;
00154   // For ease of calculation we stored interactions between types
00155   // i and j in (ni+j).  For efficiency now we coalesce the
00156   // cross interactions so that just i<=j are stored.
00157   const int arraysize = 3*pressureProfileSlabs;
00158   size_t nelems = arraysize*(numAtomTypes*(numAtomTypes+1))/2;
00159   BigReal *arr = new BigReal[nelems];
00160   memset(arr, 0, nelems*sizeof(BigReal));
00161 
00162   int i, j;
00163   for (i=0; i<numAtomTypes; i++) {
00164     for (j=0; j<numAtomTypes; j++) {
00165       int ii=i;
00166       int jj=j;
00167       if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00168       const int reductionOffset = (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00169       for (int k=0; k<arraysize; k++) {
00170         arr[reductionOffset+k] += data[k];
00171       }
00172       data += arraysize;
00173     }
00174   }
00175   // copy into reduction
00176   reduction->add(nelems, arr);
00177   delete [] arr;
00178 }
00179   
00180 void ComputeNonbondedUtil::calc_error(nonbonded *) {
00181   NAMD_bug("Tried to call missing nonbonded compute routine.");
00182 }
00183   
00184 void ComputeNonbondedUtil::select(void)
00185 {
00186   if ( CkMyRank() ) return;
00187 
00188   // These defaults die cleanly if nothing appropriate is assigned.
00189   ComputeNonbondedUtil::calcPair = calc_error;
00190   ComputeNonbondedUtil::calcPairEnergy = calc_error;
00191   ComputeNonbondedUtil::calcSelf = calc_error;
00192   ComputeNonbondedUtil::calcSelfEnergy = calc_error;
00193   ComputeNonbondedUtil::calcFullPair = calc_error;
00194   ComputeNonbondedUtil::calcFullPairEnergy = calc_error;
00195   ComputeNonbondedUtil::calcFullSelf = calc_error;
00196   ComputeNonbondedUtil::calcFullSelfEnergy = calc_error;
00197   ComputeNonbondedUtil::calcMergePair = calc_error;
00198   ComputeNonbondedUtil::calcMergePairEnergy = calc_error;
00199   ComputeNonbondedUtil::calcMergeSelf = calc_error;
00200   ComputeNonbondedUtil::calcMergeSelfEnergy = calc_error;
00201   ComputeNonbondedUtil::calcSlowPair = calc_error;
00202   ComputeNonbondedUtil::calcSlowPairEnergy = calc_error;
00203   ComputeNonbondedUtil::calcSlowSelf = calc_error;
00204   ComputeNonbondedUtil::calcSlowSelfEnergy = calc_error;
00205 
00206   SimParameters * simParams = Node::Object()->simParameters;
00207   Parameters * params = Node::Object()->parameters;
00208 
00209   table_ener = params->table_ener;
00210   rowsize = params->rowsize;
00211   columnsize = params->columnsize;
00212 
00213   commOnly = simParams->commOnly;
00214   fixedAtomsOn = ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces );
00215 
00216   cutoff = simParams->cutoff;
00217   cutoff2 = cutoff*cutoff;
00218 
00219 //fepb
00220   alchFepOn = simParams->alchFepOn;
00221   alchThermIntOn = simParams->alchThermIntOn;
00222   alchLambda = alchLambda2 = 0;
00223   lesOn = simParams->lesOn;
00224   lesScaling = lesFactor = 0;
00225   Bool tabulatedEnergies = simParams->tabulatedEnergies;
00226   alchVdwShiftCoeff = simParams->alchVdwShiftCoeff;
00227   alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
00228   alchElecLambdaStart = simParams->alchElecLambdaStart;
00229 
00230   alchDecouple = simParams->alchDecouple;
00231 
00232   delete [] lambda_table;
00233   lambda_table = 0;
00234 
00235   pairInteractionOn = simParams->pairInteractionOn;
00236   pairInteractionSelf = simParams->pairInteractionSelf;
00237   pressureProfileOn = simParams->pressureProfileOn;
00238 
00239   if ( alchFepOn ) {
00240     alchLambda = simParams->alchLambda;
00241     alchLambda2 = simParams->alchLambda2;
00242     ComputeNonbondedUtil::calcPair = calc_pair_energy_fep;
00243     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_fep;
00244     ComputeNonbondedUtil::calcSelf = calc_self_energy_fep;
00245     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_fep;
00246     ComputeNonbondedUtil::calcFullPair = calc_pair_energy_fullelect_fep;
00247     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_fep;
00248     ComputeNonbondedUtil::calcFullSelf = calc_self_energy_fullelect_fep;
00249     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_fep;
00250     ComputeNonbondedUtil::calcMergePair = calc_pair_energy_merge_fullelect_fep;
00251     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_fep;
00252     ComputeNonbondedUtil::calcMergeSelf = calc_self_energy_merge_fullelect_fep;
00253     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_fep;
00254     ComputeNonbondedUtil::calcSlowPair = calc_pair_energy_slow_fullelect_fep;
00255     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_fep;
00256     ComputeNonbondedUtil::calcSlowSelf = calc_self_energy_slow_fullelect_fep;
00257     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_fep;
00258   }  else if ( alchThermIntOn ) {
00259     alchLambda = simParams->alchLambda;
00260     ComputeNonbondedUtil::calcPair = calc_pair_ti;
00261     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_ti;
00262     ComputeNonbondedUtil::calcSelf = calc_self_ti;
00263     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_ti;
00264     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_ti;
00265     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_ti;
00266     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_ti;
00267     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_ti;
00268     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_ti;
00269     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_ti;
00270     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_ti;
00271     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_ti;
00272     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_ti;
00273     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_ti;
00274     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_ti;
00275     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_ti;
00276   } else if ( lesOn ) {
00277     lesFactor = simParams->lesFactor;
00278     lesScaling = 1.0 / (double)lesFactor;
00279     lambda_table = new BigReal[(lesFactor+1)*(lesFactor+1)];
00280     for ( int ip=0; ip<=lesFactor; ++ip ) {
00281       for ( int jp=0; jp<=lesFactor; ++jp ) {
00282         BigReal lambda_pair = 1.0;
00283         if (ip || jp ) {
00284           if (ip && jp && ip != jp) {
00285             lambda_pair = 0.0;
00286           } else {
00287             lambda_pair = lesScaling;
00288           }
00289         }
00290         lambda_table[(lesFactor+1)*ip+jp] = lambda_pair;
00291       }
00292     }
00293     ComputeNonbondedUtil::calcPair = calc_pair_les;
00294     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_les;
00295     ComputeNonbondedUtil::calcSelf = calc_self_les;
00296     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_les;
00297     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_les;
00298     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_les;
00299     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_les;
00300     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_les;
00301     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_les;
00302     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_les;
00303     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_les;
00304     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_les;
00305     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_les;
00306     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_les;
00307     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_les;
00308     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_les;
00309   } else if ( pressureProfileOn) {
00310     pressureProfileSlabs = simParams->pressureProfileSlabs;
00311     pressureProfileAtomTypes = simParams->pressureProfileAtomTypes;
00312 
00313     ComputeNonbondedUtil::calcPair = calc_pair_pprof;
00314     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_pprof;
00315     ComputeNonbondedUtil::calcSelf = calc_self_pprof;
00316     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_pprof;
00317     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_pprof;
00318     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_pprof;
00319     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_pprof;
00320     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_pprof;
00321     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_pprof;
00322     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_pprof;
00323     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_pprof;
00324     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_pprof;
00325     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_pprof;
00326     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_pprof;
00327     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_pprof;
00328     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_pprof;
00329   } else if ( pairInteractionOn ) {
00330     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_int;
00331     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_int;
00332     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_int;
00333     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_int;
00334     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_int;
00335     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_int;
00336   } else if ( tabulatedEnergies ) {
00337     ComputeNonbondedUtil::calcPair = calc_pair_tabener;
00338     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_tabener;
00339     ComputeNonbondedUtil::calcSelf = calc_self_tabener;
00340     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_tabener;
00341     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_tabener;
00342     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_tabener;
00343     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_tabener;
00344     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_tabener;
00345     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_tabener;
00346     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_tabener;
00347     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_tabener;
00348     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_tabener;
00349     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_tabener;
00350     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_tabener;
00351     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_tabener;
00352     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_tabener;
00353   } else {
00354     ComputeNonbondedUtil::calcPair = calc_pair;
00355     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy;
00356     ComputeNonbondedUtil::calcSelf = calc_self;
00357     ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy;
00358     ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect;
00359     ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect;
00360     ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect;
00361     ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect;
00362     ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect;
00363     ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect;
00364     ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect;
00365     ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect;
00366     ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect;
00367     ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect;
00368     ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect;
00369     ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect;
00370   }
00371 
00372 //fepe
00373 
00374   dielectric_1 = 1.0/simParams->dielectric;
00375   if ( ! ljTable ) ljTable = new LJTable;
00376   mol = Node::Object()->molecule;
00377   scaling = simParams->nonbondedScaling;
00378   if ( simParams->exclude == SCALED14 )
00379   {
00380     scale14 = simParams->scale14;
00381   }
00382   else
00383   {
00384     scale14 = 1.;
00385   }
00386   if ( simParams->switchingActive )
00387   {
00388     switchOn = simParams->switchingDist;
00389     switchOn_1 = 1.0/switchOn;
00390     // d0 = 1.0/(cutoff-switchOn);
00391     switchOn2 = switchOn*switchOn;
00392     c0 = 1.0/(cutoff2-switchOn2);
00393   }
00394   else
00395   {
00396     switchOn = cutoff;
00397     switchOn_1 = 1.0/switchOn;
00398     // d0 = 0.;  // avoid division by zero
00399     switchOn2 = switchOn*switchOn;
00400     c0 = 0.;  // avoid division by zero
00401   }
00402   c1 = c0*c0*c0;
00403   c3 = 3.0 * (cutoff2 - switchOn2);
00404   c5 = 0;
00405   c6 = 0;
00406   c7 = 0;
00407   c8 = 0;
00408 
00409   const int PMEOn = simParams->PMEOn;
00410 
00411   if ( PMEOn ) {
00412     ewaldcof = simParams->PMEEwaldCoefficient;
00413     BigReal TwoBySqrtPi = 1.12837916709551;
00414     pi_ewaldcof = TwoBySqrtPi * ewaldcof;
00415   }
00416 
00417   int splitType = SPLIT_NONE;
00418   if ( simParams->switchingActive ) splitType = SPLIT_SHIFT;
00419   if ( simParams->fullDirectOn || simParams->FMAOn || PMEOn ) {
00420     switch ( simParams->longSplitting ) {
00421       case C2:
00422       splitType = SPLIT_C2;
00423       break;
00424 
00425       case C1:
00426       splitType = SPLIT_C1;
00427       break;
00428 
00429       case XPLOR:
00430       NAMD_die("Sorry, XPLOR splitting not supported.");
00431       break;
00432 
00433       case SHARP:
00434       NAMD_die("Sorry, SHARP splitting not supported.");
00435       break;
00436 
00437       default:
00438       NAMD_die("Unknown splitting type found!");
00439 
00440     }
00441   }
00442 
00443   BigReal r2_tol = 0.1;
00444   
00445   r2_delta = 1.0;
00446   r2_delta_exp = 0;
00447   while ( r2_delta > r2_tol ) { r2_delta /= 2.0; r2_delta_exp += 1; }
00448   r2_delta_1 = 1.0 / r2_delta;
00449 
00450   if ( ! CkMyPe() ) {
00451     iout << iINFO << "NONBONDED TABLE R-SQUARED SPACING: " <<
00452                                 r2_delta << "\n" << endi;
00453   }
00454 
00455   BigReal r2_tmp = 1.0;
00456   int cutoff2_exp = 0;
00457   while ( (cutoff2 + r2_delta) > r2_tmp ) { r2_tmp *= 2.0; cutoff2_exp += 1; }
00458 
00459   int i;
00460   int n = (r2_delta_exp + cutoff2_exp) * 64 + 1;
00461 
00462   if ( ! CkMyPe() ) {
00463     iout << iINFO << "NONBONDED TABLE SIZE: " <<
00464                                 n << " POINTS\n" << endi;
00465   }
00466 
00467   if ( table_alloc ) delete [] table_alloc;
00468   table_alloc = new BigReal[61*n+16];
00469   BigReal *table_align = table_alloc;
00470   while ( ((long)table_align) % 128 ) ++table_align;
00471   table_noshort = table_align;
00472   table_short = table_align + 16*n;
00473   slow_table = table_align + 32*n;
00474   fast_table = table_align + 36*n;
00475   scor_table = table_align + 40*n;
00476   corr_table = table_align + 44*n;
00477   full_table = table_align + 48*n;
00478   vdwa_table = table_align + 52*n;
00479   vdwb_table = table_align + 56*n;
00480   r2_table = table_align + 60*n;
00481   BigReal *fast_i = fast_table + 4;
00482   BigReal *scor_i = scor_table + 4;
00483   BigReal *slow_i = slow_table + 4;
00484   BigReal *vdwa_i = vdwa_table + 4;
00485   BigReal *vdwb_i = vdwb_table + 4;
00486   BigReal *r2_i = r2_table;  *(r2_i++) = r2_delta;
00487   BigReal r2_limit = simParams->limitDist * simParams->limitDist;
00488   if ( r2_limit < r2_delta ) r2_limit = r2_delta;
00489   int r2_delta_i = 0;  // entry for r2 == r2_delta
00490 
00491   // fill in the table, fix up i==0 (r2==0) below
00492   for ( i=1; i<n; ++i ) {
00493 
00494     const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00495     const BigReal r2_del = r2_base / 64.0;
00496     const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00497 
00498     if ( r2 <= r2_limit ) r2_delta_i = i;
00499 
00500     const BigReal r = sqrt(r2);
00501     const BigReal r_1 = 1.0/r;
00502     const BigReal r_2 = 1.0/r2;
00503 
00504     // fast_ is defined as (full_ - slow_)
00505     // corr_ and fast_ are both zero at the cutoff, full_ is not
00506     // all three are approx 1/r at short distances
00507 
00508     // for actual interpolation, we use fast_ for fast forces and
00509     // scor_ = slow_ + corr_ - full_ and slow_ for slow forces
00510     // since these last two are of small magnitude
00511 
00512     BigReal fast_energy, fast_gradient;
00513     BigReal scor_energy, scor_gradient;
00514     BigReal slow_energy, slow_gradient;
00515 
00516     // corr_ is PME direct sum, or similar correction term
00517     // corr_energy is multiplied by r until later
00518     // corr_gradient is multiplied by -r^2 until later
00519     BigReal corr_energy, corr_gradient;
00520 
00521     
00522     if ( PMEOn ) {
00523       BigReal tmp_a = r * ewaldcof;
00524       BigReal tmp_b = erfc(tmp_a);
00525       corr_energy = tmp_b;
00526       corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
00527     } else {
00528       corr_energy = corr_gradient = 0;
00529     }
00530 
00531     switch(splitType) {
00532       case SPLIT_NONE:
00533         fast_energy = 1.0/r;
00534         fast_gradient = -1.0/r2;
00535         scor_energy = scor_gradient = 0;
00536         slow_energy = slow_gradient = 0;
00537         break;
00538       case SPLIT_SHIFT: {
00539         BigReal shiftVal = r2/cutoff2 - 1.0;
00540         shiftVal *= shiftVal;
00541         BigReal dShiftVal = 2.0 * (r2/cutoff2 - 1.0) * 2.0*r/cutoff2;
00542         fast_energy = shiftVal/r;
00543         fast_gradient = dShiftVal/r - shiftVal/r2;
00544         scor_energy = scor_gradient = 0;
00545         slow_energy = slow_gradient = 0;
00546         } 
00547         break;
00548       case SPLIT_C1:
00549         // calculate actual energy and gradient
00550         slow_energy = 0.5/cutoff * (3.0 - (r2/cutoff2));
00551         slow_gradient = -1.0/cutoff2 * (r/cutoff);
00552         // calculate scor from slow and corr
00553         scor_energy = slow_energy + (corr_energy - 1.0)/r;
00554         scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00555         // calculate fast from slow
00556         fast_energy = 1.0/r - slow_energy;
00557         fast_gradient = -1.0/r2 - slow_gradient;
00558         break;
00559       case SPLIT_C2:
00560         //
00561         // Quintic splitting function contributed by
00562         // Bruce Berne, Ruhong Zhou, and Joe Morrone
00563         //
00564         // calculate actual energy and gradient
00565         slow_energy = r2/(cutoff*cutoff2) * (6.0 * (r2/cutoff2)
00566             - 15.0*(r/cutoff) + 10.0);
00567         slow_gradient = r/(cutoff*cutoff2) * (24.0 * (r2/cutoff2)
00568             - 45.0 *(r/cutoff) + 20.0);
00569         // calculate scor from slow and corr
00570         scor_energy = slow_energy + (corr_energy - 1.0)/r;
00571         scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00572         // calculate fast from slow
00573         fast_energy = 1.0/r - slow_energy;
00574         fast_gradient = -1.0/r2 - slow_gradient;
00575         break;
00576     }
00577 
00578     // foo_gradient is calculated as ( d foo_energy / d r )
00579     // and now divided by 2r to get ( d foo_energy / d r2 )
00580 
00581     fast_gradient *= 0.5 * r_1;
00582     scor_gradient *= 0.5 * r_1;
00583     slow_gradient *= 0.5 * r_1;
00584 
00585     // let modf be 1 if excluded, 1-scale14 if modified, 0 otherwise,
00586     // add scor_ - modf * slow_ to slow terms and
00587     // add fast_ - modf * fast_ to fast terms.
00588 
00589     BigReal vdwa_energy, vdwa_gradient;
00590     BigReal vdwb_energy, vdwb_gradient;
00591 
00592     const BigReal r_6 = r_2*r_2*r_2;
00593     const BigReal r_12 = r_6*r_6;
00594 
00595     // Lennard-Jones switching function
00596     const BigReal c2 = cutoff2-r2;
00597     const BigReal c4 = c2*(c3-2.0*c2);
00598     const BigReal switchVal =         // used for Lennard-Jones
00599                         ( r2 > switchOn2 ? c2*c4*c1 : 1.0 );
00600     const BigReal dSwitchVal =        // d switchVal / d r2
00601                         ( r2 > switchOn2 ? 2*c1*(c2*c2-c4) : 0.0 );
00602 
00603     vdwa_energy = switchVal * r_12;
00604     vdwb_energy = switchVal * r_6;
00605 
00606     vdwa_gradient = ( dSwitchVal - 6.0 * switchVal * r_2 ) * r_12;
00607     vdwb_gradient = ( dSwitchVal - 3.0 * switchVal * r_2 ) * r_6;
00608 
00609 
00610     *(fast_i++) = fast_energy;
00611     *(fast_i++) = fast_gradient;
00612     *(fast_i++) = 0;
00613     *(fast_i++) = 0;
00614     *(scor_i++) = scor_energy;
00615     *(scor_i++) = scor_gradient;
00616     *(scor_i++) = 0;
00617     *(scor_i++) = 0;
00618     *(slow_i++) = slow_energy;
00619     *(slow_i++) = slow_gradient;
00620     *(slow_i++) = 0;
00621     *(slow_i++) = 0;
00622     *(vdwa_i++) = vdwa_energy;
00623     *(vdwa_i++) = vdwa_gradient;
00624     *(vdwa_i++) = 0;
00625     *(vdwa_i++) = 0;
00626     *(vdwb_i++) = vdwb_energy;
00627     *(vdwb_i++) = vdwb_gradient;
00628     *(vdwb_i++) = 0;
00629     *(vdwb_i++) = 0;
00630     *(r2_i++) = r2 + r2_delta;
00631 
00632   }
00633 
00634   if ( ! r2_delta_i ) {
00635     NAMD_bug("Failed to find table entry for r2 == r2_limit\n");
00636   }
00637   if ( r2_table[r2_delta_i] > r2_limit + r2_delta ) {
00638     NAMD_bug("Found bad table entry for r2 == r2_limit\n");
00639   }
00640 
00641   int j;
00642   const char *table_name = "XXXX";
00643   int smooth_short = 0;
00644   for ( j=0; j<5; ++j ) {
00645     BigReal *t0 = 0;
00646     switch (j) {
00647       case 0: 
00648         t0 = fast_table;
00649         table_name = "FAST";
00650         smooth_short = 1;
00651       break;
00652       case 1: 
00653         t0 = scor_table;
00654         table_name = "SCOR";
00655         smooth_short = 0;
00656       break;
00657       case 2: 
00658         t0 = slow_table;
00659         table_name = "SLOW";
00660         smooth_short = 0;
00661       break;
00662       case 3: 
00663         t0 = vdwa_table;
00664         table_name = "VDWA";
00665         smooth_short = 1;
00666       break;
00667       case 4: 
00668         t0 = vdwb_table;
00669         table_name = "VDWB";
00670         smooth_short = 1;
00671       break;
00672     }
00673     // patch up data for i=0
00674     t0[0] = t0[4] - t0[5] * ( r2_delta / 64.0 );  // energy
00675     t0[1] = t0[5];  // gradient
00676     t0[2] = 0;
00677     t0[3] = 0;
00678     if ( smooth_short ) {
00679       BigReal energy0 = t0[4*r2_delta_i];
00680       BigReal gradient0 = t0[4*r2_delta_i+1];
00681       BigReal r20 = r2_table[r2_delta_i];
00682       t0[0] = energy0 - gradient0 * (r20 - r2_table[0]);  // energy
00683       t0[1] = gradient0;  // gradient
00684     }
00685     BigReal *t;
00686     for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
00687       BigReal x = ( r2_delta * ( 1 << (i/64) ) ) / 64.0;
00688       if ( r2_table[i+1] != r2_table[i] + x ) {
00689         NAMD_bug("Bad table delta calculation.\n");
00690       }
00691       if ( smooth_short && i+1 < r2_delta_i ) {
00692         BigReal energy0 = t0[4*r2_delta_i];
00693         BigReal gradient0 = t0[4*r2_delta_i+1];
00694         BigReal r20 = r2_table[r2_delta_i];
00695         t[4] = energy0 - gradient0 * (r20 - r2_table[i+1]);  // energy
00696         t[5] = gradient0;  // gradient
00697       }
00698       BigReal v1 = t[0];
00699       BigReal g1 = t[1];
00700       BigReal v2 = t[4];
00701       BigReal g2 = t[5];
00702       // explicit formulas for v1 + g1 x + c x^2 + d x^3
00703       BigReal c = ( 3.0 * (v2 - v1) - x * (2.0 * g1 + g2) ) / ( x * x );
00704       BigReal d = ( -2.0 * (v2 - v1) + x * (g1 + g2) ) / ( x * x * x );
00705       // since v2 - v1 is imprecise, we refine c and d numerically
00706       // important because we need accurate forces (more than energies!)
00707       for ( int k=0; k < 2; ++k ) {
00708         BigReal dv = (v1 - v2) + ( ( d * x + c ) * x + g1 ) * x;
00709         BigReal dg = (g1 - g2) + ( 3.0 * d * x + 2.0 * c ) * x;
00710         c -= ( 3.0 * dv - x * dg ) / ( x * x );
00711         d -= ( -2.0 * dv + x * dg ) / ( x * x * x );
00712       }
00713       // store in the array;
00714       t[2] = c;  t[3] = d;
00715     }
00716 
00717     if ( ! CkMyPe() ) {
00718     BigReal dvmax = 0;
00719     BigReal dgmax = 0;
00720     BigReal dvmax_r = 0;
00721     BigReal dgmax_r = 0;
00722     BigReal fdvmax = 0;
00723     BigReal fdgmax = 0;
00724     BigReal fdvmax_r = 0;
00725     BigReal fdgmax_r = 0;
00726     for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
00727       const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00728       const BigReal r2_del = r2_base / 64.0;
00729       const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00730       const BigReal r = sqrt(r2);
00731       BigReal x = r2_del;
00732       BigReal dv = ( ( t[3] * x + t[2] ) * x + t[1] ) * x + t[0] - t[4];
00733       BigReal dg = ( 3.0 * t[3] * x + 2.0 * t[2] ) * x + t[1] - t[5];
00734       if ( t[4] != 0. && fabs(dv/t[4]) > fdvmax ) {
00735         fdvmax = fabs(dv/t[4]); fdvmax_r = r;
00736       }
00737       if ( fabs(dv) > dvmax ) {
00738         dvmax = fabs(dv); dvmax_r = r;
00739       }
00740       if ( t[5] != 0. && fabs(dg/t[5]) > fdgmax ) {
00741         fdgmax = fabs(dg/t[5]); fdgmax_r = r;
00742       }
00743       if ( fabs(dg) > dgmax ) {
00744         dgmax = fabs(dg); dgmax_r = r;
00745       }
00746 #if 0
00747       if (dv != 0.) CkPrintf("TABLE %d ENERGY ERROR %g AT %g (%d)\n",j,dv,r,i);
00748       if (dg != 0.) CkPrintf("TABLE %d FORCE ERROR %g AT %g (%d)\n",j,dg,r,i);
00749 #endif
00750     }
00751     if ( dvmax != 0.0 ) {
00752       iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
00753         " TABLE ENERGY: " << dvmax << " AT " << dvmax_r << "\n" << endi;
00754     }
00755     if ( fdvmax != 0.0 ) {
00756       iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
00757         " TABLE ENERGY: " << fdvmax << " AT " << fdvmax_r << "\n" << endi;
00758     }
00759     if ( dgmax != 0.0 ) {
00760       iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
00761         " TABLE FORCE: " << dgmax << " AT " << dgmax_r << "\n" << endi;
00762     }
00763     if ( fdgmax != 0.0 ) {
00764       iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
00765         " TABLE FORCE: " << fdgmax << " AT " << fdgmax_r << "\n" << endi;
00766     }
00767     }
00768 
00769   }
00770 
00771   for ( i=0; i<4*n; ++i ) {
00772     corr_table[i] = fast_table[i] + scor_table[i];
00773     full_table[i] = fast_table[i] + slow_table[i];
00774   }
00775 
00776 #if 0  
00777   for ( i=0; i<n; ++i ) {
00778    for ( int j=0; j<4; ++j ) {
00779     table_short[16*i+6-2*j] = table_noshort[16*i+6-2*j] = vdwa_table[4*i+j];
00780     table_short[16*i+7-2*j] = table_noshort[16*i+7-2*j] = vdwb_table[4*i+j];
00781     table_short[16*i+8+3-j] = fast_table[4*i+j];
00782     table_short[16*i+12+3-j] = scor_table[4*i+j];
00783     table_noshort[16*i+8+3-j] = corr_table[4*i+j];
00784     table_noshort[16*i+12+3-j] = full_table[4*i+j];
00785    }
00786   }
00787 #endif 
00788 
00789   for ( i=0; i<n; ++i ) {
00790     table_short[16*i+ 0] = table_noshort[16*i+0] = -6.*vdwa_table[4*i+3];
00791     table_short[16*i+ 2] = table_noshort[16*i+2] = -6.*vdwb_table[4*i+3];
00792     table_short[16*i+ 4] = table_noshort[16*i+4] = -2.*vdwa_table[4*i+1];
00793     table_short[16*i+ 6] = table_noshort[16*i+6] = -2.*vdwb_table[4*i+1];
00794     
00795     table_short[16*i+1] = table_noshort[16*i+1] = -4.*vdwa_table[4*i+2];
00796     table_short[16*i+3] = table_noshort[16*i+3] = -4.*vdwb_table[4*i+2];
00797     table_short[16*i+5] = table_noshort[16*i+5] = -1.*vdwa_table[4*i+0];
00798     table_short[16*i+7] = table_noshort[16*i+7] = -1.*vdwb_table[4*i+0];
00799     
00800     table_short[16*i+8]  = -6.*fast_table[4*i+3];
00801     table_short[16*i+9]  = -4.*fast_table[4*i+2];
00802     table_short[16*i+10] = -2.*fast_table[4*i+1];
00803     table_short[16*i+11] = -1.*fast_table[4*i+0];
00804 
00805     table_noshort[16*i+8]  = -6.*corr_table[4*i+3];
00806     table_noshort[16*i+9]  = -4.*corr_table[4*i+2];
00807     table_noshort[16*i+10] = -2.*corr_table[4*i+1];
00808     table_noshort[16*i+11] = -1.*corr_table[4*i+0];
00809 
00810     table_short[16*i+12] = -6.*scor_table[4*i+3];
00811     table_short[16*i+13] = -4.*scor_table[4*i+2];
00812     table_short[16*i+14] = -2.*scor_table[4*i+1];
00813     table_short[16*i+15] = -1.*scor_table[4*i+0];
00814 
00815     table_noshort[16*i+12] = -6.*full_table[4*i+3];
00816     table_noshort[16*i+13] = -4.*full_table[4*i+2];
00817     table_noshort[16*i+14] = -2.*full_table[4*i+1];
00818     table_noshort[16*i+15] = -1.*full_table[4*i+0];
00819   }
00820 
00821 #if 0
00822   char fname[100];
00823   sprintf(fname,"/tmp/namd.table.pe%d.dat",CkMyPe());
00824   FILE *f = fopen(fname,"w");
00825   for ( i=0; i<(n-1); ++i ) {
00826     const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00827     const BigReal r2_del = r2_base / 64.0;
00828     const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00829     BigReal *t;
00830     if ( r2 + r2_delta != r2_table[i] ) fprintf(f,"r2 error! ");
00831     fprintf(f,"%g",r2);
00832     t = fast_table + 4*i;
00833     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00834     t = scor_table + 4*i;
00835     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00836     t = slow_table + 4*i;
00837     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00838     t = corr_table + 4*i;
00839     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00840     t = full_table + 4*i;
00841     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00842     t = vdwa_table + 4*i;
00843     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00844     t = vdwb_table + 4*i;
00845     fprintf(f,"   %g %g %g %g", t[0], t[1], t[2], t[3]);
00846     fprintf(f,"\n");
00847   }
00848   fclose(f);
00849 #endif
00850 
00851 #ifdef NAMD_CUDA
00852   build_cuda_force_table();
00853 #endif
00854 
00855 }
00856 

Generated on Tue Nov 24 04:07:43 2009 for NAMD by  doxygen 1.3.9.1