Learn the MD of amber t3.3: implicit solvent model (GB)

One 、 brief introduction

Implicit solvent model , In this model, there are only solute molecules and no solvent molecules （ water , No ions ）, The advantage of such a model is to reduce MD Amount of computation （ After all, the number of atoms is greatly reduced ）, Moreover, the influence of solvent molecules on the potential energy change of solute can be reasonably considered .

Amber The implicit solvent model commonly used in is generalized Born Model （GB Model ）, The algorithm of this model is in the force field ff99SB and ff14Bonlysc There are good test results in .

Two 、 The process

1、 Prepare simulated proteins ：1L2Y

Download link ：https://www.rcsb.org/structure/1L2Y

For download pdb Document processing ：

pdb4amber 1l2y.pdb > xx.pdb

tleap

source leaprc.protein.ff19SB

set default PBradii mbondi3   # Guaranteed added H The atomic spacing conforms to a GB Model

protein=loadpdb xx.pdb

saveamberparm protein xx.parm7 xx.rst7

savepdb protein xx.pdb

Get simulated top and crd file , Plus H Of pdb file . 

2、 Prepare the simulated input file ： 

Minimize energy ：

We use it tleap or charmm-GUI Add water to the system , Protein plus H when , There may be bad contact between atoms . If you don't deal with this bad contact , During the simulation process, the local energy will be too high and the system will eventually be destroyed . The function of energy minimization is to find a conformation with the lowest energy , Eliminate these bad contacts .

To write 01min.in file ：

energy minimization
&cntrl
  imin = 1,                        #1 Minimize ,0 Do not minimize
  maxcyc=100,                 # Minimize at most 100 Time
  ntx = 1,                          # There has never been any speed rst Read coordinates in file
  ntwr = 100, ntpr = 10,    # The former sets the output to rst Steps , The latter setting is output to out Steps
  ioutfm=0, ntxo=1,          # The former points out nc File format binary ,netcdf, The latter designates rst File format ascii
  cut = 1000.0,                 # Non bond truncation distance （ Unit angstrom ）
  ntb=0,                            # There is no periodic boundary （ After all, there is no water ）
  igb = 8,                          # Specify that this is an implicit solvent model GBneck2

  gbsa=3, surften=0.007, # Specify solvent accessible surface area and surface tension
  saltcon = 0.0,
&end

Heat up ：

We can't directly simulate at the desired temperature , But slowly climb up from a low temperature to a high temperature ,Amber The commonly used thermostat is Langevin thermostat, Set collision constant γ(gamma=1.0ps-1), The general target temperature is 300K. In the process , We'll talk about the main chain atoms (CA,C,O,N) Add limiting force to restrain .

heating with backbone restraints
&cntrl
   ntx = 1,                            
   ntwx = 5000,ntwr = 500, ntpr = 5000,

   ntwe = 0,                  #ntwe Step output energy and temperature , Set to 0 There is no           
   ioutfm=0, ntxo=1,     # In binary format nc and ascii Format rst                
   imin = 0,
   nstlim = 500000, dt = 0.002, #dt： The length of each step , The unit is ps
   ntt = 3,gamma_ln=1.,temp0 = 100.0,   #ntt=3, use Langevin thermostat
   nscm = 1000,
   ig = -1,         # Random number seed
   ntc= 2, ntf = 2, # Yes H Constraints on the vibration of atomic bonds （shake Algorithm ）
   cut = 1000,    # Restrictions on non bonding forces between atoms , Distance greater than 1000A Just ignore the force between the two atoms

   igb=8,
   gbsa=3, surften=0.007,
   ntb = 0,
   saltcon=0.,
   ntr=1,restraintmask="@CA,N,C,O",  restraint_wt=10.0, #ntr=1： The specified atom will be constrained
   nmropt=1,    #1： Read the following &wt Area parameters                          
&end
&wt
   TYPE="TEMP0", istep1=0, istep2=500000,
   value1=100., value2=300.,        #0 Step by step , The temperature is 100K,50000 Step (1ns) The temperature is 300K    
&end
&wt
   TYPE="END",
&end 

Balance ：

In this simulation , The system is simulated at the final temperature we set , Long simulation time , There are two stages , Reduce the limiting force on the main chain atoms in turn .

To write 03equil01.in, run 250ps：

&cntrl
   ntx = 5,                              

   ntwx = 5000, ntwe = 0, ntwr = 500, ntpr = 5000,                 

   ioutfm=0, ntxo=1, 

   imin = 0, nstlim = 125000, dt = 0.002, 

   ntt = 3, gamma_ln=1.0, temp0 = 300.0,  
   nscm = 1000, ig = -1, irest = 1,        #Iirest=1： Restart the simulation , from rst Read the speed and coordinates in the file   

   ntc= 2, ntf = 2, cut = 1000,           

   igb=8, gbsa=3, surften=0.007, ntb = 0,
   saltcon=0., 

   ntr=1, restraintmask="@CA,N,C,O",     
   restraint_wt=1.0, 
/

To write 03equil02.in：

&cntrl
   ntx = 5,                             
   ntwx = 5000, ntwe = 0, ntwr = 500,    
   ntpr = 5000,                 

   ioutfm=0, ntxo=1,                     

   imin = 0, nstlim = 125000, dt = 0.002, 

   ntt = 3, gamma_ln=1., temp0 = 300.0,  
   nscm = 1000, ig = -1, irest = 1,                 

   ntc= 2, ntf = 2, cut = 1000,           

   igb=8, gbsa=3, surften=0.007, ntb = 0,
   saltcon=0., 

   ntr=1, restraintmask="@CA,N,C,O",     
   restraint_wt=0.1, 
/

  The difference between the two equilibrium simulations is restraint_wt Decrease of parameter value , from 10 To 1 To 0.1, It means that the limiting force on the main chain atoms is getting smaller and smaller , Sports will be more natural .restraint_wt Its unit is kcal/mol/angstrom2, The bigger it is , The more restricted the motion of the atom is .

3、 To simulate ：

To write sh Script ：

#！/usr/bin/bash

export CUDA_VISIBLE_DEVICES=3  # Specify to run simulated GPU model

export path=/home/packages/amber20/bin  # Description path is a global variable
$path/pmemd.cuda -O -i 01min.in -p xx.parm7 -c xx.rst7 -o min.out\
-x min.crd -inf min.info -r min.rst7

echo  -e “\033[1;5;31;40m minimize is finished \033[0m"

$path/pmemd.cuda -O -i 02heat.in -p xx.parm7 -c min.rst7 -ref min.rst7\
-o heat.out -x heat.crd -inf heat.info -r heat.rst7

echo -e “\033[1;5;31;40m heating is finished \033[0m"

$path/pmemd.cuda -O -i 03equil01.in -p xx.parm7 -c heat.rst7 -ref heat.rst7\
-o eq1.out -x eq1.crd -inf eq1.info -r eq1.rst7

echo -e “\033[1;5;31;40m equilibrium 01 is finished \033[0m"

$path/pmemd.cuda -O -i 03equil02.in -p xx.parm7 -c eq1.rst7 -ref eq1.rst7\
-o eq2.out -x eq2.crd -inf eq2.info -r eq2.rst7

echo -e “\033[1;5;31;40m equilibrium 02 is finished \033[0m"

3、 ... and 、 Reference link ：

Equilibration of Implicit Solvent