00001
00007
00008
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 Bool ComputeNonbondedUtil::commOnly;
00033 Bool ComputeNonbondedUtil::fixedAtomsOn;
00034 BigReal ComputeNonbondedUtil::cutoff;
00035 BigReal ComputeNonbondedUtil::cutoff2;
00036 BigReal ComputeNonbondedUtil::dielectric_1;
00037 const LJTable* ComputeNonbondedUtil::ljTable = 0;
00038 const Molecule* ComputeNonbondedUtil::mol;
00039 BigReal ComputeNonbondedUtil::r2_delta;
00040 BigReal ComputeNonbondedUtil::r2_delta_1;
00041 int ComputeNonbondedUtil::r2_delta_exp;
00042 int ComputeNonbondedUtil::rowsize;
00043 int ComputeNonbondedUtil::columnsize;
00044 BigReal* ComputeNonbondedUtil::table_alloc = 0;
00045 BigReal* ComputeNonbondedUtil::table_ener = 0;
00046 BigReal* ComputeNonbondedUtil::table_short;
00047 BigReal* ComputeNonbondedUtil::table_noshort;
00048 BigReal* ComputeNonbondedUtil::fast_table;
00049 BigReal* ComputeNonbondedUtil::scor_table;
00050 BigReal* ComputeNonbondedUtil::slow_table;
00051 BigReal* ComputeNonbondedUtil::corr_table;
00052 BigReal* ComputeNonbondedUtil::full_table;
00053 BigReal* ComputeNonbondedUtil::vdwa_table;
00054 BigReal* ComputeNonbondedUtil::vdwb_table;
00055 BigReal* ComputeNonbondedUtil::r2_table;
00056 BigReal ComputeNonbondedUtil::scaling;
00057 BigReal ComputeNonbondedUtil::scale14;
00058 BigReal ComputeNonbondedUtil::switchOn;
00059 BigReal ComputeNonbondedUtil::switchOn_1;
00060 BigReal ComputeNonbondedUtil::switchOn2;
00061 BigReal ComputeNonbondedUtil::v_vdwa;
00062 BigReal ComputeNonbondedUtil::v_vdwb;
00063 BigReal ComputeNonbondedUtil::k_vdwa;
00064 BigReal ComputeNonbondedUtil::k_vdwb;
00065 BigReal ComputeNonbondedUtil::cutoff_3;
00066 BigReal ComputeNonbondedUtil::cutoff_6;
00067 BigReal ComputeNonbondedUtil::c0;
00068 BigReal ComputeNonbondedUtil::c1;
00069 BigReal ComputeNonbondedUtil::c3;
00070 BigReal ComputeNonbondedUtil::c5;
00071 BigReal ComputeNonbondedUtil::c6;
00072 BigReal ComputeNonbondedUtil::c7;
00073 BigReal ComputeNonbondedUtil::c8;
00074
00075
00076 Bool ComputeNonbondedUtil::alchFepOn;
00077 Bool ComputeNonbondedUtil::alchThermIntOn;
00078 Bool ComputeNonbondedUtil::Fep_WCA_repuOn;
00079 Bool ComputeNonbondedUtil::Fep_WCA_dispOn;
00080 BigReal ComputeNonbondedUtil::WCA_rcut1;
00081 BigReal ComputeNonbondedUtil::WCA_rcut2;
00082 BigReal ComputeNonbondedUtil::alchLambda;
00083 BigReal ComputeNonbondedUtil::alchLambda2;
00084 BigReal ComputeNonbondedUtil::alchVdwShiftCoeff;
00085 BigReal ComputeNonbondedUtil::alchVdwLambdaEnd;
00086 BigReal ComputeNonbondedUtil::alchElecLambdaStart;
00087 Bool ComputeNonbondedUtil::alchDecouple;
00088
00089 Bool ComputeNonbondedUtil::lesOn;
00090 int ComputeNonbondedUtil::lesFactor;
00091 BigReal ComputeNonbondedUtil::lesScaling;
00092
00093 BigReal* ComputeNonbondedUtil::lambda_table = 0;
00094
00095 Bool ComputeNonbondedUtil::pairInteractionOn;
00096 Bool ComputeNonbondedUtil::pairInteractionSelf;
00097
00098 Bool ComputeNonbondedUtil::pressureProfileOn;
00099 int ComputeNonbondedUtil::pressureProfileSlabs;
00100 int ComputeNonbondedUtil::pressureProfileAtomTypes;
00101 BigReal ComputeNonbondedUtil::pressureProfileThickness;
00102 BigReal ComputeNonbondedUtil::pressureProfileMin;
00103
00104 Bool ComputeNonbondedUtil::accelMDOn;
00105
00106 Bool ComputeNonbondedUtil::drudeNbthole;
00107
00108 BigReal ComputeNonbondedUtil::ewaldcof;
00109 BigReal ComputeNonbondedUtil::pi_ewaldcof;
00110
00111
00112 Bool ComputeNonbondedUtil::goForcesOn;
00113 int ComputeNonbondedUtil::goMethod;
00114
00115
00116 void (*ComputeNonbondedUtil::calcPair)(nonbonded *);
00117 void (*ComputeNonbondedUtil::calcPairEnergy)(nonbonded *);
00118 void (*ComputeNonbondedUtil::calcSelf)(nonbonded *);
00119 void (*ComputeNonbondedUtil::calcSelfEnergy)(nonbonded *);
00120
00121 void (*ComputeNonbondedUtil::calcFullPair)(nonbonded *);
00122 void (*ComputeNonbondedUtil::calcFullPairEnergy)(nonbonded *);
00123 void (*ComputeNonbondedUtil::calcFullSelf)(nonbonded *);
00124 void (*ComputeNonbondedUtil::calcFullSelfEnergy)(nonbonded *);
00125
00126 void (*ComputeNonbondedUtil::calcMergePair)(nonbonded *);
00127 void (*ComputeNonbondedUtil::calcMergePairEnergy)(nonbonded *);
00128 void (*ComputeNonbondedUtil::calcMergeSelf)(nonbonded *);
00129 void (*ComputeNonbondedUtil::calcMergeSelfEnergy)(nonbonded *);
00130
00131 void (*ComputeNonbondedUtil::calcSlowPair)(nonbonded *);
00132 void (*ComputeNonbondedUtil::calcSlowPairEnergy)(nonbonded *);
00133 void (*ComputeNonbondedUtil::calcSlowSelf)(nonbonded *);
00134 void (*ComputeNonbondedUtil::calcSlowSelfEnergy)(nonbonded *);
00135
00136
00137 #define SPLIT_NONE 1
00138 #define SPLIT_SHIFT 2
00139 #define SPLIT_C1 3
00140 #define SPLIT_XPLOR 4
00141 #define SPLIT_C2 5
00142 #define SPLIT_MARTINI 6
00143
00144 void ComputeNonbondedUtil::submitReductionData(BigReal *data, SubmitReduction *reduction)
00145 {
00146 reduction->item(REDUCTION_EXCLUSION_CHECKSUM) += data[exclChecksumIndex];
00147 reduction->item(REDUCTION_PAIRLIST_WARNINGS) += data[pairlistWarningIndex];
00148 reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
00149 reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += data[fullElectEnergyIndex];
00150 reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
00151
00152 reduction->item(REDUCTION_GO_NATIVE_ENERGY) += data[goNativeEnergyIndex];
00153 reduction->item(REDUCTION_GO_NONNATIVE_ENERGY) += data[goNonnativeEnergyIndex];
00154
00155
00156 reduction->item(REDUCTION_ELECT_ENERGY_F) += data[electEnergyIndex_s];
00157 reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += data[fullElectEnergyIndex_s];
00158 reduction->item(REDUCTION_LJ_ENERGY_F) += data[vdwEnergyIndex_s];
00159
00160 reduction->item(REDUCTION_ELECT_ENERGY_TI_1) += data[electEnergyIndex_ti_1];
00161 reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1) += data[fullElectEnergyIndex_ti_1];
00162 reduction->item(REDUCTION_LJ_ENERGY_TI_1) += data[vdwEnergyIndex_ti_1];
00163 reduction->item(REDUCTION_ELECT_ENERGY_TI_2) += data[electEnergyIndex_ti_2];
00164 reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2) += data[fullElectEnergyIndex_ti_2];
00165 reduction->item(REDUCTION_LJ_ENERGY_TI_2) += data[vdwEnergyIndex_ti_2];
00166
00167 ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00168 ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,fullElectVirialIndex);
00169 ADD_VECTOR(reduction,REDUCTION_PAIR_VDW_FORCE,data,pairVDWForceIndex);
00170 ADD_VECTOR(reduction,REDUCTION_PAIR_ELECT_FORCE,data,pairElectForceIndex);
00171 reduction->item(REDUCTION_COMPUTE_CHECKSUM) += 1.;
00172 }
00173
00174 void ComputeNonbondedUtil::submitPressureProfileData(BigReal *data,
00175 SubmitReduction *reduction)
00176 {
00177 if (!reduction) return;
00178 int numAtomTypes = pressureProfileAtomTypes;
00179
00180
00181
00182 const int arraysize = 3*pressureProfileSlabs;
00183 size_t nelems = arraysize*(numAtomTypes*(numAtomTypes+1))/2;
00184 BigReal *arr = new BigReal[nelems];
00185 memset(arr, 0, nelems*sizeof(BigReal));
00186
00187 int i, j;
00188 for (i=0; i<numAtomTypes; i++) {
00189 for (j=0; j<numAtomTypes; j++) {
00190 int ii=i;
00191 int jj=j;
00192 if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00193 const int reductionOffset = (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00194 for (int k=0; k<arraysize; k++) {
00195 arr[reductionOffset+k] += data[k];
00196 }
00197 data += arraysize;
00198 }
00199 }
00200
00201 reduction->add(nelems, arr);
00202 delete [] arr;
00203 }
00204
00205 void ComputeNonbondedUtil::calc_error(nonbonded *) {
00206 NAMD_bug("Tried to call missing nonbonded compute routine.");
00207 }
00208
00209 void ComputeNonbondedUtil::select(void)
00210 {
00211 if ( CkMyRank() ) return;
00212
00213
00214 ComputeNonbondedUtil::calcPair = calc_error;
00215 ComputeNonbondedUtil::calcPairEnergy = calc_error;
00216 ComputeNonbondedUtil::calcSelf = calc_error;
00217 ComputeNonbondedUtil::calcSelfEnergy = calc_error;
00218 ComputeNonbondedUtil::calcFullPair = calc_error;
00219 ComputeNonbondedUtil::calcFullPairEnergy = calc_error;
00220 ComputeNonbondedUtil::calcFullSelf = calc_error;
00221 ComputeNonbondedUtil::calcFullSelfEnergy = calc_error;
00222 ComputeNonbondedUtil::calcMergePair = calc_error;
00223 ComputeNonbondedUtil::calcMergePairEnergy = calc_error;
00224 ComputeNonbondedUtil::calcMergeSelf = calc_error;
00225 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_error;
00226 ComputeNonbondedUtil::calcSlowPair = calc_error;
00227 ComputeNonbondedUtil::calcSlowPairEnergy = calc_error;
00228 ComputeNonbondedUtil::calcSlowSelf = calc_error;
00229 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_error;
00230
00231 SimParameters * simParams = Node::Object()->simParameters;
00232 Parameters * params = Node::Object()->parameters;
00233
00234 table_ener = params->table_ener;
00235 rowsize = params->rowsize;
00236 columnsize = params->columnsize;
00237
00238 commOnly = simParams->commOnly;
00239 fixedAtomsOn = ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces );
00240
00241 cutoff = simParams->cutoff;
00242 cutoff2 = cutoff*cutoff;
00243
00244
00245 alchFepOn = simParams->alchFepOn;
00246 Fep_WCA_repuOn = simParams->alchFepWCARepuOn;
00247 Fep_WCA_dispOn = simParams->alchFepWCADispOn;
00248 alchThermIntOn = simParams->alchThermIntOn;
00249 alchLambda = alchLambda2 = 0;
00250 lesOn = simParams->lesOn;
00251 lesScaling = lesFactor = 0;
00252 Bool tabulatedEnergies = simParams->tabulatedEnergies;
00253 alchVdwShiftCoeff = simParams->alchVdwShiftCoeff;
00254 WCA_rcut1 = simParams->alchFepWCArcut1;
00255 WCA_rcut2 = simParams->alchFepWCArcut2;
00256 alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
00257 alchElecLambdaStart = simParams->alchElecLambdaStart;
00258
00259 alchDecouple = simParams->alchDecouple;
00260
00261 delete [] lambda_table;
00262 lambda_table = 0;
00263
00264 pairInteractionOn = simParams->pairInteractionOn;
00265 pairInteractionSelf = simParams->pairInteractionSelf;
00266 pressureProfileOn = simParams->pressureProfileOn;
00267
00268
00269 goForcesOn = simParams->goForcesOn;
00270 goMethod = simParams->goMethod;
00271
00272
00273 accelMDOn = simParams->accelMDOn;
00274
00275 drudeNbthole = simParams->drudeOn && (simParams->drudeNbtholeCut > 0.0);
00276
00277 if ( drudeNbthole ) {
00278 #ifdef NAMD_CUDA
00279 NAMD_die("drudeNbthole is not supported in CUDA version");
00280 #endif
00281 if ( alchFepOn )
00282 NAMD_die("drudeNbthole is not supported with alchemical free-energy perturbation");
00283 if ( alchThermIntOn )
00284 NAMD_die("drudeNbthole is not supported with alchemical thermodynamic integration");
00285 if ( lesOn )
00286 NAMD_die("drudeNbthole is not supported with locally enhanced sampling");
00287 if ( pairInteractionOn )
00288 NAMD_die("drudeNbthole is not supported with pair interaction calculation");
00289 if ( pressureProfileOn )
00290 NAMD_die("drudeNbthole is not supported with pressure profile calculation");
00291 }
00292
00293 if ( alchFepOn ) {
00294 #ifdef NAMD_CUDA
00295 NAMD_die("Alchemical free-energy perturbation is not supported in CUDA version");
00296 #endif
00297 alchLambda = simParams->alchLambda;
00298 alchLambda2 = simParams->alchLambda2;
00299 ComputeNonbondedUtil::calcPair = calc_pair_energy_fep;
00300 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_fep;
00301 ComputeNonbondedUtil::calcSelf = calc_self_energy_fep;
00302 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_fep;
00303 ComputeNonbondedUtil::calcFullPair = calc_pair_energy_fullelect_fep;
00304 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_fep;
00305 ComputeNonbondedUtil::calcFullSelf = calc_self_energy_fullelect_fep;
00306 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_fep;
00307 ComputeNonbondedUtil::calcMergePair = calc_pair_energy_merge_fullelect_fep;
00308 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_fep;
00309 ComputeNonbondedUtil::calcMergeSelf = calc_self_energy_merge_fullelect_fep;
00310 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_fep;
00311 ComputeNonbondedUtil::calcSlowPair = calc_pair_energy_slow_fullelect_fep;
00312 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_fep;
00313 ComputeNonbondedUtil::calcSlowSelf = calc_self_energy_slow_fullelect_fep;
00314 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_fep;
00315 } else if ( alchThermIntOn ) {
00316 #ifdef NAMD_CUDA
00317 NAMD_die("Alchemical thermodynamic integration is not supported in CUDA version");
00318 #endif
00319 alchLambda = simParams->alchLambda;
00320 ComputeNonbondedUtil::calcPair = calc_pair_ti;
00321 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_ti;
00322 ComputeNonbondedUtil::calcSelf = calc_self_ti;
00323 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_ti;
00324 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_ti;
00325 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_ti;
00326 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_ti;
00327 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_ti;
00328 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_ti;
00329 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_ti;
00330 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_ti;
00331 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_ti;
00332 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_ti;
00333 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_ti;
00334 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_ti;
00335 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_ti;
00336 } else if ( lesOn ) {
00337 #ifdef NAMD_CUDA
00338 NAMD_die("Locally enhanced sampling is not supported in CUDA version");
00339 #endif
00340 lesFactor = simParams->lesFactor;
00341 lesScaling = 1.0 / (double)lesFactor;
00342 lambda_table = new BigReal[(lesFactor+1)*(lesFactor+1)];
00343 for ( int ip=0; ip<=lesFactor; ++ip ) {
00344 for ( int jp=0; jp<=lesFactor; ++jp ) {
00345 BigReal lambda_pair = 1.0;
00346 if (ip || jp ) {
00347 if (ip && jp && ip != jp) {
00348 lambda_pair = 0.0;
00349 } else {
00350 lambda_pair = lesScaling;
00351 }
00352 }
00353 lambda_table[(lesFactor+1)*ip+jp] = lambda_pair;
00354 }
00355 }
00356 ComputeNonbondedUtil::calcPair = calc_pair_les;
00357 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_les;
00358 ComputeNonbondedUtil::calcSelf = calc_self_les;
00359 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_les;
00360 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_les;
00361 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_les;
00362 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_les;
00363 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_les;
00364 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_les;
00365 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_les;
00366 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_les;
00367 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_les;
00368 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_les;
00369 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_les;
00370 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_les;
00371 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_les;
00372 } else if ( pressureProfileOn) {
00373 #ifdef NAMD_CUDA
00374 NAMD_die("Pressure profile calculation is not supported in CUDA version");
00375 #endif
00376 pressureProfileSlabs = simParams->pressureProfileSlabs;
00377 pressureProfileAtomTypes = simParams->pressureProfileAtomTypes;
00378
00379 ComputeNonbondedUtil::calcPair = calc_pair_pprof;
00380 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_pprof;
00381 ComputeNonbondedUtil::calcSelf = calc_self_pprof;
00382 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_pprof;
00383 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_pprof;
00384 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_pprof;
00385 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_pprof;
00386 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_pprof;
00387 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_pprof;
00388 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_pprof;
00389 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_pprof;
00390 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_pprof;
00391 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_pprof;
00392 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_pprof;
00393 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_pprof;
00394 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_pprof;
00395 } else if ( pairInteractionOn ) {
00396 #ifdef NAMD_CUDA
00397 NAMD_die("Pair interaction calculation is not supported in CUDA version");
00398 #endif
00399 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_int;
00400 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_int;
00401 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_int;
00402 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_int;
00403 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_int;
00404 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_int;
00405 } else if ( tabulatedEnergies ) {
00406 #ifdef NAMD_CUDA
00407 NAMD_die("Tabulated energies is not supported in CUDA version");
00408 #endif
00409 ComputeNonbondedUtil::calcPair = calc_pair_tabener;
00410 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_tabener;
00411 ComputeNonbondedUtil::calcSelf = calc_self_tabener;
00412 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_tabener;
00413 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_tabener;
00414 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_tabener;
00415 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_tabener;
00416 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_tabener;
00417 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_tabener;
00418 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_tabener;
00419 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_tabener;
00420 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_tabener;
00421 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_tabener;
00422 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_tabener;
00423 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_tabener;
00424 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_tabener;
00425 } else if ( goForcesOn ) {
00426 #ifdef NAMD_CUDA
00427 NAMD_die("Go forces is not supported in CUDA version");
00428 #endif
00429 ComputeNonbondedUtil::calcPair = calc_pair_go;
00430 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_go;
00431 ComputeNonbondedUtil::calcSelf = calc_self_go;
00432 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy_go;
00433 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect_go;
00434 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect_go;
00435 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect_go;
00436 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect_go;
00437 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect_go;
00438 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect_go;
00439 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect_go;
00440 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect_go;
00441 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect_go;
00442 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect_go;
00443 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect_go;
00444 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect_go;
00445 } else {
00446 ComputeNonbondedUtil::calcPair = calc_pair;
00447 ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy;
00448 ComputeNonbondedUtil::calcSelf = calc_self;
00449 ComputeNonbondedUtil::calcSelfEnergy = calc_self_energy;
00450 ComputeNonbondedUtil::calcFullPair = calc_pair_fullelect;
00451 ComputeNonbondedUtil::calcFullPairEnergy = calc_pair_energy_fullelect;
00452 ComputeNonbondedUtil::calcFullSelf = calc_self_fullelect;
00453 ComputeNonbondedUtil::calcFullSelfEnergy = calc_self_energy_fullelect;
00454 ComputeNonbondedUtil::calcMergePair = calc_pair_merge_fullelect;
00455 ComputeNonbondedUtil::calcMergePairEnergy = calc_pair_energy_merge_fullelect;
00456 ComputeNonbondedUtil::calcMergeSelf = calc_self_merge_fullelect;
00457 ComputeNonbondedUtil::calcMergeSelfEnergy = calc_self_energy_merge_fullelect;
00458 ComputeNonbondedUtil::calcSlowPair = calc_pair_slow_fullelect;
00459 ComputeNonbondedUtil::calcSlowPairEnergy = calc_pair_energy_slow_fullelect;
00460 ComputeNonbondedUtil::calcSlowSelf = calc_self_slow_fullelect;
00461 ComputeNonbondedUtil::calcSlowSelfEnergy = calc_self_energy_slow_fullelect;
00462 }
00463
00464
00465
00466 dielectric_1 = 1.0/simParams->dielectric;
00467 if ( ! ljTable ) ljTable = new LJTable;
00468 mol = Node::Object()->molecule;
00469 scaling = simParams->nonbondedScaling;
00470 if ( simParams->exclude == SCALED14 )
00471 {
00472 scale14 = simParams->scale14;
00473 }
00474 else
00475 {
00476 scale14 = 1.;
00477 }
00478 if ( simParams->switchingActive )
00479 {
00480 switchOn = simParams->switchingDist;
00481 switchOn_1 = 1.0/switchOn;
00482
00483 switchOn2 = switchOn*switchOn;
00484 c0 = 1.0/(cutoff2-switchOn2);
00485
00486 if ( simParams->vdwForceSwitching ) {
00487 double switchOn3 = switchOn * switchOn2;
00488 double cutoff3 = cutoff * cutoff2;
00489 double switchOn6 = switchOn3 * switchOn3;
00490 double cutoff6 = cutoff3 * cutoff3;
00491 v_vdwa = -1. / ( switchOn6 * cutoff6 );
00492 v_vdwb = -1. / ( switchOn3 * cutoff3 );
00493 k_vdwa = cutoff6 / ( cutoff6 - switchOn6 );
00494 k_vdwb = cutoff3 / ( cutoff3 - switchOn3 );
00495 cutoff_3 = 1. / cutoff3;
00496 cutoff_6 = 1. / cutoff6;
00497 }
00498 }
00499 else
00500 {
00501 switchOn = cutoff;
00502 switchOn_1 = 1.0/switchOn;
00503
00504 switchOn2 = switchOn*switchOn;
00505 c0 = 0.;
00506 }
00507 c1 = c0*c0*c0;
00508 c3 = 3.0 * (cutoff2 - switchOn2);
00509 c5 = 0;
00510 c6 = 0;
00511 c7 = 0;
00512 c8 = 0;
00513
00514 const int PMEOn = simParams->PMEOn;
00515 const int MSMOn = simParams->MSMOn;
00516 const int MSMSplit = simParams->MSMSplit;
00517
00518 if ( PMEOn ) {
00519 ewaldcof = simParams->PMEEwaldCoefficient;
00520 BigReal TwoBySqrtPi = 1.12837916709551;
00521 pi_ewaldcof = TwoBySqrtPi * ewaldcof;
00522 }
00523
00524 int splitType = SPLIT_NONE;
00525 if ( simParams->switchingActive ) splitType = SPLIT_SHIFT;
00526 if ( simParams->martiniSwitching ) splitType = SPLIT_MARTINI;
00527 if ( simParams->fullDirectOn || simParams->FMAOn || PMEOn || MSMOn ) {
00528 switch ( simParams->longSplitting ) {
00529 case C2:
00530 splitType = SPLIT_C2;
00531 break;
00532
00533 case C1:
00534 splitType = SPLIT_C1;
00535 break;
00536
00537 case XPLOR:
00538 NAMD_die("Sorry, XPLOR splitting not supported.");
00539 break;
00540
00541 case SHARP:
00542 NAMD_die("Sorry, SHARP splitting not supported.");
00543 break;
00544
00545 default:
00546 NAMD_die("Unknown splitting type found!");
00547
00548 }
00549 }
00550
00551 BigReal r2_tol = 0.1;
00552
00553 r2_delta = 1.0;
00554 r2_delta_exp = 0;
00555 while ( r2_delta > r2_tol ) { r2_delta /= 2.0; r2_delta_exp += 1; }
00556 r2_delta_1 = 1.0 / r2_delta;
00557
00558 if ( ! CkMyPe() ) {
00559 iout << iINFO << "NONBONDED TABLE R-SQUARED SPACING: " <<
00560 r2_delta << "\n" << endi;
00561 }
00562
00563 BigReal r2_tmp = 1.0;
00564 int cutoff2_exp = 0;
00565 while ( (cutoff2 + r2_delta) > r2_tmp ) { r2_tmp *= 2.0; cutoff2_exp += 1; }
00566
00567 int i;
00568 int n = (r2_delta_exp + cutoff2_exp) * 64 + 1;
00569
00570 if ( ! CkMyPe() ) {
00571 iout << iINFO << "NONBONDED TABLE SIZE: " <<
00572 n << " POINTS\n" << endi;
00573 }
00574
00575 if ( table_alloc ) delete [] table_alloc;
00576 table_alloc = new BigReal[61*n+16];
00577 BigReal *table_align = table_alloc;
00578 while ( ((long)table_align) % 128 ) ++table_align;
00579 table_noshort = table_align;
00580 table_short = table_align + 16*n;
00581 slow_table = table_align + 32*n;
00582 fast_table = table_align + 36*n;
00583 scor_table = table_align + 40*n;
00584 corr_table = table_align + 44*n;
00585 full_table = table_align + 48*n;
00586 vdwa_table = table_align + 52*n;
00587 vdwb_table = table_align + 56*n;
00588 r2_table = table_align + 60*n;
00589 BigReal *fast_i = fast_table + 4;
00590 BigReal *scor_i = scor_table + 4;
00591 BigReal *slow_i = slow_table + 4;
00592 BigReal *vdwa_i = vdwa_table + 4;
00593 BigReal *vdwb_i = vdwb_table + 4;
00594 BigReal *r2_i = r2_table; *(r2_i++) = r2_delta;
00595 BigReal r2_limit = simParams->limitDist * simParams->limitDist;
00596 if ( r2_limit < r2_delta ) r2_limit = r2_delta;
00597 int r2_delta_i = 0;
00598
00599
00600 for ( i=1; i<n; ++i ) {
00601
00602 const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00603 const BigReal r2_del = r2_base / 64.0;
00604 const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00605
00606 if ( r2 <= r2_limit ) r2_delta_i = i;
00607
00608 const BigReal r = sqrt(r2);
00609 const BigReal r_1 = 1.0/r;
00610 const BigReal r_2 = 1.0/r2;
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 BigReal fast_energy, fast_gradient;
00621 BigReal scor_energy, scor_gradient;
00622 BigReal slow_energy, slow_gradient;
00623
00624
00625
00626
00627 BigReal corr_energy, corr_gradient;
00628
00629
00630 if ( PMEOn ) {
00631 BigReal tmp_a = r * ewaldcof;
00632 BigReal tmp_b = erfc(tmp_a);
00633 corr_energy = tmp_b;
00634 corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
00635 } else if ( MSMOn ) {
00636 BigReal a_1 = 1.0/cutoff;
00637 BigReal r_a = r * a_1;
00638 BigReal g, dg;
00639 SPOLY(&g, &dg, r_a, MSMSplit);
00640 corr_energy = 1 - r_a * g;
00641 corr_gradient = 1 + r_a*r_a * dg;
00642 } else {
00643 corr_energy = corr_gradient = 0;
00644 }
00645
00646 switch(splitType) {
00647 case SPLIT_NONE:
00648 fast_energy = 1.0/r;
00649 fast_gradient = -1.0/r2;
00650 scor_energy = scor_gradient = 0;
00651 slow_energy = slow_gradient = 0;
00652 break;
00653 case SPLIT_SHIFT: {
00654 BigReal shiftVal = r2/cutoff2 - 1.0;
00655 shiftVal *= shiftVal;
00656 BigReal dShiftVal = 2.0 * (r2/cutoff2 - 1.0) * 2.0*r/cutoff2;
00657 fast_energy = shiftVal/r;
00658 fast_gradient = dShiftVal/r - shiftVal/r2;
00659 scor_energy = scor_gradient = 0;
00660 slow_energy = slow_gradient = 0;
00661 }
00662 break;
00663 case SPLIT_MARTINI: {
00664
00665 const BigReal COUL_SWITCH = 0.;
00666
00667 const BigReal p1 = 1.;
00668 BigReal A1 = p1 * ((p1+1)*COUL_SWITCH-(p1+4)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,2));
00669 BigReal B1 = -p1 * ((p1+1)*COUL_SWITCH-(p1+3)*cutoff)/(pow(cutoff,p1+2)*pow(cutoff-COUL_SWITCH,3));
00670 BigReal X1 = 1.0/pow(cutoff,p1)-A1/3.0*pow(cutoff-COUL_SWITCH,3)-B1/4.0*pow(cutoff-COUL_SWITCH,4);
00671 BigReal r12 = (r-COUL_SWITCH)*(r-COUL_SWITCH);
00672 BigReal r13 = (r-COUL_SWITCH)*(r-COUL_SWITCH)*(r-COUL_SWITCH);
00673 BigReal shiftVal = -(A1/3.0)*r13 - (B1/4.0)*r12*r12 - X1;
00674 BigReal dShiftVal = A1*r12 + B1*r13;
00675 fast_energy = (1/r) + shiftVal;
00676 fast_gradient = 1/(r2) + dShiftVal;
00677 scor_energy = scor_gradient = 0;
00678 slow_energy = slow_gradient = 0;
00679 }
00680 break;
00681 case SPLIT_C1:
00682
00683 slow_energy = 0.5/cutoff * (3.0 - (r2/cutoff2));
00684 slow_gradient = -1.0/cutoff2 * (r/cutoff);
00685
00686 scor_energy = slow_energy + (corr_energy - 1.0)/r;
00687 scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00688
00689 fast_energy = 1.0/r - slow_energy;
00690 fast_gradient = -1.0/r2 - slow_gradient;
00691 break;
00692 case SPLIT_C2:
00693
00694
00695
00696
00697
00698 slow_energy = r2/(cutoff*cutoff2) * (6.0 * (r2/cutoff2)
00699 - 15.0*(r/cutoff) + 10.0);
00700 slow_gradient = r/(cutoff*cutoff2) * (24.0 * (r2/cutoff2)
00701 - 45.0 *(r/cutoff) + 20.0);
00702
00703 scor_energy = slow_energy + (corr_energy - 1.0)/r;
00704 scor_gradient = slow_gradient - (corr_gradient - 1.0)/r2;
00705
00706 fast_energy = 1.0/r - slow_energy;
00707 fast_gradient = -1.0/r2 - slow_gradient;
00708 break;
00709 }
00710
00711
00712
00713
00714 fast_gradient *= 0.5 * r_1;
00715 scor_gradient *= 0.5 * r_1;
00716 slow_gradient *= 0.5 * r_1;
00717
00718
00719
00720
00721
00722 BigReal vdwa_energy, vdwa_gradient;
00723 BigReal vdwb_energy, vdwb_gradient;
00724
00725 const BigReal r_6 = r_2*r_2*r_2;
00726 const BigReal r_12 = r_6*r_6;
00727
00728
00729 if ( simParams->vdwForceSwitching ) {
00730
00731 if ( r2 > switchOn2 ) {
00732 BigReal tmpa = r_6 - cutoff_6;
00733 vdwa_energy = k_vdwa * tmpa * tmpa;
00734 BigReal tmpb = r_1 * r_2 - cutoff_3;
00735 vdwb_energy = k_vdwb * tmpb * tmpb;
00736 vdwa_gradient = -6.0 * k_vdwa * tmpa * r_2 * r_6;
00737 vdwb_gradient = -3.0 * k_vdwb * tmpb * r_2 * r_2 * r_1;
00738 } else {
00739 vdwa_energy = r_12 + v_vdwa;
00740 vdwb_energy = r_6 + v_vdwb;
00741 vdwa_gradient = -6.0 * r_2 * r_12;
00742 vdwb_gradient = -3.0 * r_2 * r_6;
00743 }
00744 } else if ( simParams->martiniSwitching ) {
00745
00746 BigReal r12 = (r-switchOn)*(r-switchOn); BigReal r13 = (r-switchOn)*(r-switchOn)*(r-switchOn);
00747
00748 BigReal p6 = 6;
00749 BigReal A6temp = p6 * ((p6+1)*switchOn-(p6+4)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,2));
00750 BigReal B6temp = -p6 * ((p6+1)*switchOn-(p6+3)*cutoff)/(pow(cutoff,p6+2)*pow(cutoff-switchOn,3));
00751 BigReal C6temp = 1.0/pow(cutoff,p6)-A6temp/3.0*pow(cutoff-switchOn,3)-B6temp/4.0*pow(cutoff-switchOn,4); BigReal A6 = -A6temp;
00752 BigReal B6 = -B6temp; BigReal C6 = -C6temp;
00753
00754 BigReal p12 = 12;
00755 BigReal A12 = p12 * ((p12+1)*switchOn-(p12+4)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,2));
00756 BigReal B12 = -p12 * ((p12+1)*switchOn-(p12+3)*cutoff)/(pow(cutoff,p12+2)*pow(cutoff-switchOn,3));
00757 BigReal C12 = 1.0/pow(cutoff,p12)-A12/3.0*pow(cutoff-switchOn,3)-B12/4.0*pow(cutoff-switchOn,4);
00758
00759 BigReal LJshifttempA = -(A12/3)*r13 - (B12/4)*r12*r12 - C12;
00760 BigReal LJshifttempB = -(A6/3)*r13 - (B6/4)*r12*r12 - C6;
00761 const BigReal shiftValA =
00762 ( r2 > switchOn2 ? LJshifttempA : -C12);
00763 const BigReal shiftValB =
00764 ( r2 > switchOn2 ? LJshifttempB : -C6);
00765
00766 BigReal LJdshifttempA = A12*r12 + B12*r13;
00767 BigReal LJdshifttempB = A6*r12 + B6*r13;
00768 const BigReal dshiftValA =
00769 ( r2 > switchOn2 ? LJdshifttempA*0.5*r_1 : 0 );
00770 const BigReal dshiftValB =
00771 ( r2 > switchOn2 ? LJdshifttempB*0.5*r_1 : 0 );
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781 vdwa_energy = r_12 - shiftValA;
00782 vdwb_energy = r_6 - shiftValB;
00783
00784 vdwa_gradient = -6/pow(r,14) + dshiftValA ;
00785 vdwb_gradient = -3/pow(r,8) + dshiftValB;
00786
00787 } else {
00788 const BigReal c2 = cutoff2-r2;
00789 const BigReal c4 = c2*(c3-2.0*c2);
00790 const BigReal switchVal =
00791 ( r2 > switchOn2 ? c2*c4*c1 : 1.0 );
00792 const BigReal dSwitchVal =
00793 ( r2 > switchOn2 ? 2*c1*(c2*c2-c4) : 0.0 );
00794
00795 vdwa_energy = switchVal * r_12;
00796 vdwb_energy = switchVal * r_6;
00797
00798 vdwa_gradient = ( dSwitchVal - 6.0 * switchVal * r_2 ) * r_12;
00799 vdwb_gradient = ( dSwitchVal - 3.0 * switchVal * r_2 ) * r_6;
00800 }
00801
00802
00803 *(fast_i++) = fast_energy;
00804 *(fast_i++) = fast_gradient;
00805 *(fast_i++) = 0;
00806 *(fast_i++) = 0;
00807 *(scor_i++) = scor_energy;
00808 *(scor_i++) = scor_gradient;
00809 *(scor_i++) = 0;
00810 *(scor_i++) = 0;
00811 *(slow_i++) = slow_energy;
00812 *(slow_i++) = slow_gradient;
00813 *(slow_i++) = 0;
00814 *(slow_i++) = 0;
00815 *(vdwa_i++) = vdwa_energy;
00816 *(vdwa_i++) = vdwa_gradient;
00817 *(vdwa_i++) = 0;
00818 *(vdwa_i++) = 0;
00819 *(vdwb_i++) = vdwb_energy;
00820 *(vdwb_i++) = vdwb_gradient;
00821 *(vdwb_i++) = 0;
00822 *(vdwb_i++) = 0;
00823 *(r2_i++) = r2 + r2_delta;
00824
00825 }
00826
00827 if ( ! r2_delta_i ) {
00828 NAMD_bug("Failed to find table entry for r2 == r2_limit\n");
00829 }
00830 if ( r2_table[r2_delta_i] > r2_limit + r2_delta ) {
00831 NAMD_bug("Found bad table entry for r2 == r2_limit\n");
00832 }
00833
00834 int j;
00835 const char *table_name = "XXXX";
00836 int smooth_short = 0;
00837 for ( j=0; j<5; ++j ) {
00838 BigReal *t0 = 0;
00839 switch (j) {
00840 case 0:
00841 t0 = fast_table;
00842 table_name = "FAST";
00843 smooth_short = 1;
00844 break;
00845 case 1:
00846 t0 = scor_table;
00847 table_name = "SCOR";
00848 smooth_short = 0;
00849 break;
00850 case 2:
00851 t0 = slow_table;
00852 table_name = "SLOW";
00853 smooth_short = 0;
00854 break;
00855 case 3:
00856 t0 = vdwa_table;
00857 table_name = "VDWA";
00858 smooth_short = 1;
00859 break;
00860 case 4:
00861 t0 = vdwb_table;
00862 table_name = "VDWB";
00863 smooth_short = 1;
00864 break;
00865 }
00866
00867 t0[0] = t0[4] - t0[5] * ( r2_delta / 64.0 );
00868 t0[1] = t0[5];
00869 t0[2] = 0;
00870 t0[3] = 0;
00871 if ( smooth_short ) {
00872 BigReal energy0 = t0[4*r2_delta_i];
00873 BigReal gradient0 = t0[4*r2_delta_i+1];
00874 BigReal r20 = r2_table[r2_delta_i];
00875 t0[0] = energy0 - gradient0 * (r20 - r2_table[0]);
00876 t0[1] = gradient0;
00877 }
00878 BigReal *t;
00879 for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
00880 BigReal x = ( r2_delta * ( 1 << (i/64) ) ) / 64.0;
00881 if ( r2_table[i+1] != r2_table[i] + x ) {
00882 NAMD_bug("Bad table delta calculation.\n");
00883 }
00884 if ( smooth_short && i+1 < r2_delta_i ) {
00885 BigReal energy0 = t0[4*r2_delta_i];
00886 BigReal gradient0 = t0[4*r2_delta_i+1];
00887 BigReal r20 = r2_table[r2_delta_i];
00888 t[4] = energy0 - gradient0 * (r20 - r2_table[i+1]);
00889 t[5] = gradient0;
00890 }
00891 BigReal v1 = t[0];
00892 BigReal g1 = t[1];
00893 BigReal v2 = t[4];
00894 BigReal g2 = t[5];
00895
00896 BigReal c = ( 3.0 * (v2 - v1) - x * (2.0 * g1 + g2) ) / ( x * x );
00897 BigReal d = ( -2.0 * (v2 - v1) + x * (g1 + g2) ) / ( x * x * x );
00898
00899
00900 for ( int k=0; k < 2; ++k ) {
00901 BigReal dv = (v1 - v2) + ( ( d * x + c ) * x + g1 ) * x;
00902 BigReal dg = (g1 - g2) + ( 3.0 * d * x + 2.0 * c ) * x;
00903 c -= ( 3.0 * dv - x * dg ) / ( x * x );
00904 d -= ( -2.0 * dv + x * dg ) / ( x * x * x );
00905 }
00906
00907 t[2] = c; t[3] = d;
00908 }
00909
00910 if ( ! CkMyPe() ) {
00911 BigReal dvmax = 0;
00912 BigReal dgmax = 0;
00913 BigReal dvmax_r = 0;
00914 BigReal dgmax_r = 0;
00915 BigReal fdvmax = 0;
00916 BigReal fdgmax = 0;
00917 BigReal fdvmax_r = 0;
00918 BigReal fdgmax_r = 0;
00919 for ( i=0,t=t0; i<(n-1); ++i,t+=4 ) {
00920 const BigReal r2_base = r2_delta * ( 1 << (i/64) );
00921 const BigReal r2_del = r2_base / 64.0;
00922 const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
00923 const BigReal r = sqrt(r2);
00924 BigReal x = r2_del;
00925 BigReal dv = ( ( t[3] * x + t[2] ) * x + t[1] ) * x + t[0] - t[4];
00926 BigReal dg = ( 3.0 * t[3] * x + 2.0 * t[2] ) * x + t[1] - t[5];
00927 if ( t[4] != 0. && fabs(dv/t[4]) > fdvmax ) {
00928 fdvmax = fabs(dv/t[4]); fdvmax_r = r;
00929 }
00930 if ( fabs(dv) > dvmax ) {
00931 dvmax = fabs(dv); dvmax_r = r;
00932 }
00933 if ( t[5] != 0. && fabs(dg/t[5]) > fdgmax ) {
00934 fdgmax = fabs(dg/t[5]); fdgmax_r = r;
00935 }
00936 if ( fabs(dg) > dgmax ) {
00937 dgmax = fabs(dg); dgmax_r = r;
00938 }
00939 #if 0
00940 if (dv != 0.) CkPrintf("TABLE %d ENERGY ERROR %g AT %g (%d)\n",j,dv,r,i);
00941 if (dg != 0.) CkPrintf("TABLE %d FORCE ERROR %g AT %g (%d)\n",j,dg,r,i);
00942 #endif
00943 }
00944 if ( dvmax != 0.0 ) {
00945 iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
00946 " TABLE ENERGY: " << dvmax << " AT " << dvmax_r << "\n" << endi;
00947 }
00948 if ( fdvmax != 0.0 ) {
00949 iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
00950 " TABLE ENERGY: " << fdvmax << " AT " << fdvmax_r << "\n" << endi;
00951 }
00952 if ( dgmax != 0.0 ) {
00953 iout << iINFO << "ABSOLUTE IMPRECISION IN " << table_name <<
00954 " TABLE FORCE: " << dgmax << " AT " << dgmax_r << "\n" << endi;
00955 }
00956 if ( fdgmax != 0.0 ) {
00957 iout << iINFO << "RELATIVE IMPRECISION IN " << table_name <<
00958 " TABLE FORCE: " << fdgmax << " AT " << fdgmax_r << "\n" << endi;
00959 }
00960 }
00961
00962 }
00963
00964 for ( i=0; i<4*n; ++i ) {
00965 corr_table[i] = fast_table[i] + scor_table[i];
00966 full_table[i] = fast_table[i] + slow_table[i];
00967 }
00968
00969 #if 0
00970 for ( i=0; i<n; ++i ) {
00971 for ( int j=0; j<4; ++j ) {
00972 table_short[16*i+6-2*j] = table_noshort[16*i+6-2*j] = vdwa_table[4*i+j];
00973 table_short[16*i+7-2*j] = table_noshort[16*i+7-2*j] = vdwb_table[4*i+j];
00974 table_short[16*i+8+3-j] = fast_table[4*i+j];
00975 table_short[16*i+12+3-j] = scor_table[4*i+j];
00976 table_noshort[16*i+8+3-j] = corr_table[4*i+j];
00977 table_noshort[16*i+12+3-j] = full_table[4*i+j];
00978 }
00979 }
00980 #endif
00981
00982 for ( i=0; i<n; ++i ) {
00983 table_short[16*i+ 0] = table_noshort[16*i+0] = -6.*vdwa_table[4*i+3];
00984 table_short[16*i+ 2] = table_noshort[16*i+2] = -6.*vdwb_table[4*i+3];
00985 table_short[16*i+ 4] = table_noshort[16*i+4] = -2.*vdwa_table[4*i+1];
00986 table_short[16*i+ 6] = table_noshort[16*i+6] = -2.*vdwb_table[4*i+1];
00987
00988 table_short[16*i+1] = table_noshort[16*i+1] = -4.*vdwa_table[4*i+2];
00989 table_short[16*i+3] = table_noshort[16*i+3] = -4.*vdwb_table[4*i+2];
00990 table_short[16*i+5] = table_noshort[16*i+5] = -1.*vdwa_table[4*i+0];
00991 table_short[16*i+7] = table_noshort[16*i+7] = -1.*vdwb_table[4*i+0];
00992
00993 table_short[16*i+8] = -6.*fast_table[4*i+3];
00994 table_short[16*i+9] = -4.*fast_table[4*i+2];
00995 table_short[16*i+10] = -2.*fast_table[4*i+1];
00996 table_short[16*i+11] = -1.*fast_table[4*i+0];
00997
00998 table_noshort[16*i+8] = -6.*corr_table[4*i+3];
00999 table_noshort[16*i+9] = -4.*corr_table[4*i+2];
01000 table_noshort[16*i+10] = -2.*corr_table[4*i+1];
01001 table_noshort[16*i+11] = -1.*corr_table[4*i+0];
01002
01003 table_short[16*i+12] = -6.*scor_table[4*i+3];
01004 table_short[16*i+13] = -4.*scor_table[4*i+2];
01005 table_short[16*i+14] = -2.*scor_table[4*i+1];
01006 table_short[16*i+15] = -1.*scor_table[4*i+0];
01007
01008 table_noshort[16*i+12] = -6.*full_table[4*i+3];
01009 table_noshort[16*i+13] = -4.*full_table[4*i+2];
01010 table_noshort[16*i+14] = -2.*full_table[4*i+1];
01011 table_noshort[16*i+15] = -1.*full_table[4*i+0];
01012 }
01013
01014 #if 0
01015 char fname[100];
01016 sprintf(fname,"/tmp/namd.table.pe%d.dat",CkMyPe());
01017 FILE *f = fopen(fname,"w");
01018 for ( i=0; i<(n-1); ++i ) {
01019 const BigReal r2_base = r2_delta * ( 1 << (i/64) );
01020 const BigReal r2_del = r2_base / 64.0;
01021 const BigReal r2 = r2_base - r2_delta + r2_del * (i%64);
01022 BigReal *t;
01023 if ( r2 + r2_delta != r2_table[i] ) fprintf(f,"r2 error! ");
01024 fprintf(f,"%g",r2);
01025 t = fast_table + 4*i;
01026 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01027 t = scor_table + 4*i;
01028 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01029 t = slow_table + 4*i;
01030 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01031 t = corr_table + 4*i;
01032 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01033 t = full_table + 4*i;
01034 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01035 t = vdwa_table + 4*i;
01036 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01037 t = vdwb_table + 4*i;
01038 fprintf(f," %g %g %g %g", t[0], t[1], t[2], t[3]);
01039 fprintf(f,"\n");
01040 }
01041 fclose(f);
01042 #endif
01043
01044 #ifdef NAMD_CUDA
01045 send_build_cuda_force_table();
01046 #endif
01047
01048 }
01049