Update/fix incomplete behavior of nonbondedScaling 72/4172/1
authorradakb <brian.radak@gmail.com>
Mon, 14 May 2018 20:27:20 +0000 (16:27 -0400)
committerradakb <brian.radak@gmail.com>
Mon, 14 May 2018 20:27:20 +0000 (16:27 -0400)
nonbondedScaling is supposed to permit application of a positive
scalar multiple to all ELECT and VDW terms. Since LJcorrection terms
are part of VDW, these should be also be scaled (also the virial
correction, in analogy to the force).

This is a real edge case "bug," I only realized it when trying to
perform simulated scaling on a box of rigid water, which should be
identical to simulated tempering. Even then the results only differ
appreciably when the effective temperature is >320 K. Needless to
say I'm probably the only person to ever attempt this.

TODO:

There are probably also missing ELECT terms when using Drude, the
Thole terms in particular. Since those are implemented as a 4-body
bonded term, this should be implemented easily using the "scale"
attribute of the tuple, only I don't know where that is set. Help?

If there are no objections, I would also like to add a
dihedralScaling keyword that serves the same function for (proper)
dihedrals and crossterms (maybe only when mergeCrossTerms is True?)

We could also add an option for a step-dependent update for both of
these keywords, just as we do for alchLamdba via alchLambda2 and
alchLambdaFreq. That would enable simple enhanced sampling
integrators based on nonequilibrium switches. This should offer less
of a challenge than implementing solute scaling switches at the C
level.

Change-Id: I5ea29616da95a517fee56c69a5daa87ed4edbdb7

src/Molecule.C

index 79739a4..3c77000 100644 (file)
@@ -10082,25 +10082,28 @@ void Molecule::compute_LJcorrection() {
 
 // Convenience function to simplify lambda scaling.
 BigReal Molecule::getEnergyTailCorr(const BigReal alchLambda){
+  const BigReal scl = simParams->nonbondedScaling;
   if (simParams->alchOn && simParams->alchVdwLambdaEnd) {
     const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
     const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
-    return vdw_lambda_1*tail_corr_dUdl_1 + vdw_lambda_2*tail_corr_dUdl_2;
+    return scl*(vdw_lambda_1*tail_corr_dUdl_1 + vdw_lambda_2*tail_corr_dUdl_2);
   }
   else {
-    return tail_corr_ener;
+    return scl*tail_corr_ener;
   }
 }
 
 // Convenience function to simplify lambda scaling.
 BigReal Molecule::getVirialTailCorr(const BigReal alchLambda){
+  const BigReal scl = simParams->nonbondedScaling;
   if (simParams->alchOn && simParams->alchVdwLambdaEnd) {
     const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
     const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
-    return vdw_lambda_1*tail_corr_virial_1 + vdw_lambda_2*tail_corr_virial_2;
+    return scl*(vdw_lambda_1*tail_corr_virial_1 +
+                vdw_lambda_2*tail_corr_virial_2);
   }
   else {
-    return tail_corr_virial;
+    return scl*tail_corr_virial;
   }
 }
 #endif