 Home Theory Analysis Examples GaMD-Amber GaMD-NAMD GaMD-GENESIS Contact
 Manual of GaMD in Amber A PDF Manual of latest GaMD, LiGaMD and Pep-GaMD in Amber: GaMD_Amber-manual.pdf 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)        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 ◇