Lambda-ABF implementation for the NAMD interface

[jhenin]

As fully documented in the Lagardère et al. paper
complements the Tinker-HP implementation, which is maintained on the Tinker-HP repo.

TODO

• Regression tests
• Document current limitations in ref. manual

Update: why are there no regression tests? Because of a synchronization issue in NAMD, where the current alchemical colvar interface is plugged into a "non-urgent" output routine, which causes it to run in a variable order with respect to the Colvars calculation. This causes variations in the instantaneous trajectory, even though it cancels out of the ensemble averages, so the free energies are correct. This causes failures in regression tests, which are by design sensitive to trajectory details. The GPU-resident code path of NAMD is immune from this problem and produces reproducible trajectories: it is therefore the recommended way to use this code in NAMD.

[fhh2626]

Hi @jhenin , I just read your excellent λ-ABF paper and then found this PR. I can't wait to give it a try. Is there any guide about how to write the Colvars file?
Thanks,
Haohao

[jhenin]

Hi Haohao,
The NAMD implementation still needs some refinement: it comes with performance limitations, namely, you cannot use MTS and you need energy computation every time step. Other than that, it works perfectly, and the performance is not so bad at all.

Here is an example input (using HMR):

timestep                2.0
fullElectFrequency      1
nonbondedFreq           1

# For NAMD3
computeEnergies 1

alch                    on
alchFile                        $alchconfig
alchCol                 B                    
alchOutFreq             1000

alchOutFile             $outName.fepout
alchElecLambdaStart     0.5                  
alchVdwLambdaEnd        1.0            
alchVdwShiftCoeff       5.0         
alchdecouple            on

alchType                TI
alchlambda 0

colvars                on 

cv config "
smp off  # Necessary for shared ABF
colvar {
  name l
  extendedLagrangian on
  extendedMass 150000
  extendedLangevinDamping 1000

   lowerBoundary 0
   upperBoundary 1
   reflectingLowerBoundary
   reflectingUpperBoundary
   width 0.01
  alchLambda {
  }
  outputTotalForce
  outputAppliedforce
}

 abf {
   colvars l
   fullSamples 1000
   outputFreq  50000
   historyFreq 100000
   shared on
   sharedFreq 5000
}
"

cv configfile $CVconfig
cv configfile common/protein_tilt.colvars

cv config "
colvarsTrajFrequency 100
colvarsRestartFrequency 100000"

[fhh2626]

Hi Haohao, The NAMD implementation still needs some refinement: it comes with performance limitations, namely, you cannot use MTS and you need energy computation every time step. Other than that, it works perfectly, and the performance is not so bad at all.

Here is an example input (using HMR):

timestep                2.0
fullElectFrequency      1
nonbondedFreq           1

# For NAMD3
computeEnergies 1

alch                    on
alchFile                        $alchconfig
alchCol                 B                    
alchOutFreq             1000

alchOutFile             $outName.fepout
alchElecLambdaStart     0.5                  
alchVdwLambdaEnd        1.0            
alchVdwShiftCoeff       5.0         
alchdecouple            on

alchType                TI
alchlambda 0

colvars                on 

cv config "
smp off  # Necessary for shared ABF
colvar {
  name l
  extendedLagrangian on
  extendedMass 150000
  extendedLangevinDamping 1000

   lowerBoundary 0
   upperBoundary 1
   reflectingLowerBoundary
   reflectingUpperBoundary
   width 0.01
  alchLambda {
  }
  outputTotalForce
  outputAppliedforce
}

 abf {
   colvars l
   fullSamples 1000
   outputFreq  50000
   historyFreq 100000
   shared on
   sharedFreq 5000
}
"

cv configfile $CVconfig
cv configfile common/protein_tilt.colvars

cv config "
colvarsTrajFrequency 100
colvarsRestartFrequency 100000"

Thank you very much! Just one additional question. I need to use HMR to enable a timestep of 2fs?

[jhenin]

In alchemical simulations, I usually scale down the timestep, because one of the end points has a molecule in gas phase. Depending on the nature of the ligand, you may get stable simulations at larger timesteps.

[fhh2626]

Thanks for your clarification!

[fhh2626]

Hi, @jhenin ,
I suggest also make available the Colvars harmonic restraints as the alchemical CVs, which will let us completely give up FEP/TI in binding free energy calculations.
Haohao

[jhenin]

Hi, @jhenin , I suggest also make available the Colvars harmonic restraints as the alchemical CVs, which will let us completely give up FEP/TI in binding free energy calculations. Haohao

Thanks Haohao, I like the idea.

[giacomofiorin]

Looks good overall save for some minor comments.

[giacomofiorin]

I would be okay with this in a pinch, but ideally I'd like to distinguish it from other "external" objects that are computed in a reproducible way (centers of mass, volumetric maps). Are you using this just to circumvent the check below?

        cvm::error("Error: the calculated value of colvar \""+name+
                   "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
                   "last read from the state file:\n"+cvm::to_str(x_restart)+
                   "\nPossible causes are changes in configuration, "
                   "wrong state file, or how PBC wrapping is handled.\n",
                   COLVARS_INPUT_ERROR);

[jhenin]

It is used in many places in colvar.cpp. Another way to describe this feature is: this colvar is not a function of atom coordinates, it contains a special CVC that is communicated to and from the back-end through a separate mechanism. I'm happy to rename it if you have a better name, but I don't see an alternative to this feature.

[giacomofiorin]

What do you mean by "driven alchemical parameters"?

[jhenin]

"Driven" means Colvars is responsible for integrating (driving) the dynamics of this parameter. For the current version of this, see

colvars/src/colvar.cpp

Line 332 in 5e37f65

if (!is_enabled(f_cv_total_force_current_step)) {