Manual of GaMD in Amber

GaMD Algorithm

GaMD scheme

Credit: Clarisse Ricci

GaMD {
     If (irest_gamd == 0) then
        For i = 1, …, ntcmd // run initial conventional molecular dynamics
           If (i >= ntcmdprep) Update Vmax, Vmin
           If (i >= ntcmdprep && i%ntave ==0) Update Vavg, sigmaV
        End
       Save Vmax,Vmin,Vavg,sigmaV to “gamd_restart.dat” file
       Calc_E_k0(iE,sigma0,Vmax,Vmin,Vavg,sigmaV)

       For i = ntcmd+1, …, ntcmd+nteb // Run biasing molecular dynamics simulation steps
             deltaV = 0.5*k0*(E-V)**2/(Vmax-Vmin)
             V = V + deltaV
             If (i >= ntcmd+ntebprep) Update Vmax, Vmin
             If (i >= ntcmd+ntebprep && i%ntave ==0) Update Vavg, sigmaV
             Calc_E_k0(iE,sigma0,Vmax,Vmin,Vavg,sigmaV)
       End
       Save Vmax,Vmin,Vavg,sigmaV to “gamd_restart.dat” file
     else if (irest_gamd == 1) then
       Read Vmax,Vmin,Vavg,sigmaV from “gamd_restart.dat” file
     End if
    
     For i = ntcmd+nteb+1, …, nstlim // run production simulation
       deltaV = 0.5*k0*(E-V)**2/(Vmax-Vmin)
       V = V + deltaV
     End
}

Subroutine Calc_E_k0(iE,sigma0,Vmax,Vmin,Vavg,sigmaV) {
if iE = 1 :
         E = Vmax
         k0’ = (sigma0/sigmaV) * (Vmax-Vmin)/(Vmax-Vavg)
         k0 = min(1.0, k0’)
else if iE = 2 :
         k0” = (1-sigma0/sigmaV) * (Vmax-Vmin)/(Vavg-Vmin)
         if 0 < k0” <= 1 :
                        k0 = k0”
                        E = Vmin + (Vmax-Vmin)/k0
         else
                        E = Vmax
                        k0’ = (sigma0/sigmaV) * (Vmax-Vmin)/(Vmax-Vavg)
                        k0 = min(1.0, k0’)
         end
end
}

Simulation Parameter Settings

* igamd = 0 (default)
igamd = 0: no boost is applied
igamd = 1: boost on the total potential energy only
igamd = 2: boost on the dihedral energy only
igamd = 3: boost on both dihedral and total potential energy

* iE = 1 (default)
iE = 1: set the threshold energy to the lower bound E = Vmax
iE = 2: set the threshold energy to the upper bound E = Vmin + (Vmax - Vmin)/k0

* ntcmdprep = 200,000 (default)
The number of preparation conventional molecular dynamics steps. This is used for system equilibration and the potential energies are not collected for calculating their statistics. The default is 200,000 for a simulation with 2 fs timestep.

* ntcmd = 1000000 (default)
The number of initial conventional molecular dynamics simulation steps used to
calculate the maximum, minimum, average and standard deviation of the system
potential energies (i.e., Vmax, Vmin, Vavg, sigmaV). The default is 1,000,000 for a
simulation with 2 fs timestep.

* ntebprep = 200,000 (default)
The number of preparation biasing molecular dynamics simulation steps. This is used for system equilibration after adding the boost potential and the potential statistics (Vmax, Vmin, Vavg, V) are not updated during these steps. The default is 200,000 for a simulation with 2 fs timestep.

* nteb = 1000000 (default)
The number of biasing molecular dynamics simulation steps. Potential statistics (Vmax, Vmin, Vavg, V) are updated between the ntebprep and nteb steps and used to calculate the GaMD acceleration parameters, particularly E and k0. The default is 1,000,000 for a simulation with 2 fs timestep. A greater value may be needed to ensure that the potential statistics and GaMD acceleration parameters level off before production simulation.

* ntave = 50000 (default)
The number of simulation steps used to calculate the average and standard deviation of potential energies. This variable has already been used in Amber. The default is set to 50,000 for GaMD simulations. It is recommended to be updated as about 4 times of the total number of atoms in the system. Note that ntcmd and nteb need to be multiples of ntave. The ntcmd and nteb simulation frames are recommended to be excluded for final trajectory analysis.

* sigma0P = 6.0 ! 10*kB*T
The upper limit of the standard deviation of the total potential boost that allows
for accurate reweighting if igamd is set to 1 or 3. The default is 6.0 (unit: kcal/mol).

* sigma0D = 6.0 ! 10*kB*T
The upper limit of the standard deviation of the dihedral potential boost that
allows for accurate reweighting if igamd is set to 2 or 3. The default is 6.0 (unit:
kcal/mol).
* irest_gamd = 0 (default)
irest_gamd = 0: new GaMD simulation. A file "gamd-restart.dat" that stores the potential
Vmax, Vmin, Vavg, sigmaV values will be saved automatically for simulation restart.
irest_gamd = 1: restart GaMD simulation. The "gamd-restart.dat" file will be read for restart.

GaMD-Amber Manual:

A PDF Manual of GaMD in Amber: GaMD_Amber-manual.pdf

Known Problems & Solutions

* For new GaMD simulation, there appears abrupt change in the dihedral and total potential energies when "irest=1" (reading both coordinates and velocities) is used. This leads to unexpected large boost potential (last two columns in gamd.log file), for which the normal values should be below 50 kcal/mol.

Solution: Before the bug is fixed, in the first job please try to use “irest=0, ntx=1” (reading only cooridates). After that, feel free to use irest = 1 for simulation restart.

 

Copyright Reserved © Yinglong Miao
Last updated: Tue, November 22, 2016