GaMD-Amber: Manual

Download a PDF Manual of the latest GaMD, LiGaMD, Pep-GaMD and PPI-GaMD in Amber.

GaMD Algorithm

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)

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+nteb+1, …, nstlim // run production simulation
deltaV =
0.5*k0*(E-V)**2/(Vmax-Vmin)
V = V + deltaV
End
}

PPI-GaMD Algorithm

PPI-GaMD {
If (irest_gamd == 0) then
For i = 1, …, ntcmd // run initial conventional molecular dynamics
If (i >= ntcmdprep) Update Vmax and Vmin of interaction potential energy
If (i >= ntcmdprep && i%ntave ==0) Update Vavg and sigmaV of interaction potential energy
End
Save Vmax,Vmin,Vavg,sigmaV of interaction potential energy 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 and Vmin of interaction potential energy
If (i >= ntcmd+ntebprep && i%ntave ==0) Update Vavg and sigmaV of interaction potential energy
Calc_E_k0(iE,sigma0,Vmax,Vmin,Vavg,sigmaV)
End
Save Vmax,Vmin,Vavg and sigmaV of of interaction potential energy to “gamd_restart.dat” file
else if (irest_gamd == 1) then
Read Vmax,Vmin,Vavg and sigmaV of interaction potential energy 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.

bgpro2atm Start atomic number of the second protein.

edpro2atm The final atomic number of the second protein.

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.