Manual of GaMD in Amber

A PDF Manual of latest GaMD, LiGaMD and Pep-GaMD in Amber: GaMD_Amber-manual.pdf

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
}

LiGaMD Algorithm

LiGaMD {
     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 
     
     ntwin = ntmsd+nftau*ntwx
     lig0=1 // ID of the bound ligand

     For i = ntcmd+nteb+1, …, nstlim // run production simulation
       If (ibblig == 1 && i%ntwx ==0) then // identify the bound ligand according to the distance to protein
         For ilig = 1, …, nlig
            dlig = distance(atom_p, atom_l)
            If (dmin <= dlig) blig_min=ilig; dmin=dlig
         End
         If (dmin <= dblig) blig=blig_min
       else if (ibblig == 2 && i%ntwin ==0) then // identify the bound ligand according to MSD
         For ilig = 1, …, nlig
            dlig = msd(atom_l, ntmsd, nftau)
            If (dmin <= dlig) blig_min=ilig; dmin=dlig
         End
         If (dmin <= dblig) blig=blig_min
       End if
       If (blig != lig0) Swap atomic coordinates, forces and velocities of ligand blig with lig0 for selective higher boost
     
       deltaV = 0.5*k0*(E-V)**2/(Vmax-Vmin)
       V = V + deltaV
     End
}

Pep-GaMD Algorithm

Pep-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
}

Simulation Parameter Settings

igamd Flag to apply boost potential
= 0 (default) no boost is applied
= 1 boost on the total potential energy only
= 2 boost on the dihedral energy only
= 3 dual boost on both dihedral and total potential energy
= 4 boost on the non-bonded potential energy only
= 5 dual boost on both dihedral and non-bonded potential energy
= 10 boost on non-bonded potential energy of selected region (defined by timask1 and scmask1) as for a ligand (LiGaMD)
= 11 dual boost on both non-bonded potential energy of the bound ligand and remaining potential energy of the rest of the system (LiGaMD_Dual)
= 14 boost on the total potential energy of selected region (defined by timask1 and scmask1) as for a peptide (Pep-GaMD)
= 15 dual boost on both the peptide potential energy and the total system potential energy other than the peptide potential energy (Pep-GaMD_Dual)


iE Flag to set the threshold energy E for applying all boost potentials
= 1 (default) set the threshold energy to the lower bound E = Vmax
= 2 set the threshold energy to the upper bound E = Vmin + (Vmax - Vmin)/k0


iEP Flag to overwrite iE and set the threshold energy E for applying the first boost potential in dual-boost schemes
= 1 (default) set the threshold energy to the lower bound E = Vmax
= 2 set the threshold energy to the upper bound E = Vmin + (Vmax - Vmin)/k0


iED Flag to overwrite iE and set the threshold energy E for applying the second boost potential in dual-boost schemes
= 1 (default) set the threshold energy to the lower bound E = Vmax
= 2 set the threshold energy to the upper bound E = Vmin + (Vmax - Vmin)/k0


ntcmdprep 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 The number of initial conventional molecular dynamics simulation steps. Potential energies are collected between ntcmdprep and ntcmd to calculate their maximum, minimum, average and standard deviation (Vmax, Vmin, Vavg, sigmaV). The default is 1,000,000 for a simulation with 2 fs timestep.


ntebprep 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, sigmaV) are not updated during these steps. The default is 200,000 for a simulation with 2 fs timestep.


nteb The number of biasing molecular dynamics simulation steps. Potential statistics (Vmax, Vmin, Vavg, sigmaV) 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 running production simulation between the nteb and nstlim (total simulation length) steps. Moreover, nteb can be set to nstlim, by which the potential statistics and GaMD acceleration parameters are updated adaptively throughout the simulation. This in some cases provides more appropriate acceleration.


ntave 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.


irest_gamd Flag to restart GaMD simulation
= 0 (default) new simulation. A file "gamd-restart.dat" that stores the maximum, minimum, average and standard deviation of the potential energies needed to calculate the boost potentials (depending on the igamd flag) will be saved automatically after GaMD equilibration stage.
= 1 restart simulation (ntcmd and nteb are set to 0 in this case). The "gamd-restart.dat" file will be read for restart.


sigma0P The upper limit of the standard deviation of the first potential boost that allows for accurate reweighting. The default is 6.0 (unit: kcal/mol).


sigma0D The upper limit of the standard deviation of the second potential boost that allows for accurate reweighting in dual-boost simulations (e.g., igamd = 2, 3, 5, 11 and 15). The default is 6.0 (unit: kcal/mol).


timask1 Specifies atoms of the first (bound) ligand or peptide in ambmask format when igamd = 10, 11, 14 or 15. The default is an empty string.


scmask1 Specifies atoms of the first (bound) ligand that will be described using soft core in ambmask format in LiGaMD when igamd = 10 or 11. In Pep-GaMD with igamd = 14 or 15, this flag was only used to specify atoms of peptide in ambmask format, but the peptide atoms will not be described using soft core. The default is an empty string.


nlig The total number of ligand molecules in the system. The default is 0. For LiGaMD simulations with nlig>1, it is important to note that: (1) To prepare the simulation starting structure, put the bound ligand first in the PDB, followed by the unbound ligand molecules. (2) Use residue ID of the first (bound) ligand in timask1 and scmask1, for which higher boost potential will be selectively applied. More flexible settings will be provided in future developments of the code.


ibblig The flag to boost the bound ligand selectively with nlig > 1
= 0 (default) no selective boost
= 1 boost the bound ligand selectively out of nlig ligand molecules in the system
= 2 boost the bound ligand selectively out of nlig ligand molecules in the system based on mean-square displacements (MSD)


dblig (Optional) The cutoff distance between atoms atom_p and atom_l used to identify the ligand that is bound in the protein when ibblig = 1 or the cutoff MSD of atom atom_l used to identify the ligand that is bound in the protein when ibblig = 2. If dblig (default 1.0d99 angstroms) is not set in the input file, the first boost potential will be selectively applied to the ligand with the smallest distance to the protein (ibblig = 1) or the smallest MSD (ibblig = 2) out of nlig ligand molecules in the system.


atom_p Serial number of a protein atom (starting from 1 for the first protein atom) used to calculate the ligand distance. It is used only when ibblig = 1. The default is 0.


atom_l Serial number of a ligand atom (starting from 1 for the first ligand atom) used to calculate the ligand distance to the protein. It is used only when ibblig = 1 or 2. The default is 0.


ntmsd Number of timesteps corresponding to the lagging time used to calculate the ligand MSD. It is used only when ibblig = 2. The default is 50000.


nftau Number of saved simulation frames used to calculate the ligand MSD. MSD is calculated for every time window of ntwin = ntmsd + nftau*ntwx steps, for which simulation frames are saved every ntwx steps. It is used only when ibblig = 2. The default is 10.

 

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: This bug has been fixed in AMBER 20. 

 

Copyright Reserved © Miao Lab, University of Kansas
◇ Last updated: Mon, July 18, 2022