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

• [x] Regression tests

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.