Creating Input Files:
edit tleap.in

loadAmberParams frcmod.ionsjc_tip3p(if this is not loaded, there will be error "Could not find vdW (or other) parameters for type: Cl-" ) source leaprc.constph # Load the PDBs 6pti = loadPDB 6pti_fixed.pdb # Make the disulfide bonds bond 6pti.5.SG 6pti.55.SG bond 6pti.14.SG 6pti.38.SG bond 6pti.30.SG 6pti.51.SG # Solvate solvateBox 6pti TIP3PBOX 10.0 #Add ions(make the salt concentration 0.1M) addIonsRand 6pti CL 13 NA 7 # Save topology files saveAmberParm 6pti 6pti.parm7 6pti.rst7 # Quit quit

tleap -f tleap.in

Creating Cpin File

cpinutil.py -p 6pti.parm7 -o 6pti.cpin -op 6pti.mod0.parm7

note: edit line 256 in cpinutil.py

changeradii---->changeRadii

-op: Carboxylate residues in explicit solvent simulations require a modified topology file, thus need to print a new parameter file.

System Minimization

Minimization to prepare for constant pH MD &cntrl imin = 1,                                   ! Turn on minimization ncyc=1000,                              ! If ncyc<maxcyc steepest descent before switching to conjugated gradient for remaining maxcyc = 5000,                       ! Total number of minimization cycles cut = 8,                                    ! Nonbonded cutoff in angstroms ntpr=50,                                   ! Print frequency restraint_wt = 10,                     ! 10 kcal/mol/A**2 restraint force constant restraintmask = '@CA,C,O,N'     ! Restraints on the backbone atoms only /

mpirun -np 12 pmemd.MPI -O -i min.mdin -p 6pti.mod0.parm7 -c 6pti.rst7 -o 6pti.min.mdout -r 6pti.min.rst7 -ref 6pti.rst7 -cpin 6pti.cpin

Heating The System

Explicit solvent constant pH initial heating mdin &cntrl imin=0, irest=0, ntx=1, ntpr=1000, ntwx=1000, nstlim=200000, dt=0.002, ntt=3, gamma_ln=5.0, ig=-1, ntc=2,ntf=2,cut=8,ntb=1, iwrap=1,ioutfm=1,nmropt=1, / &wt TYPE='TEMP0', ISTEP1=0, ISTEP2=150000, VALUE1=10.0, VALUE2=300.0, / &wt TYPE='END' /

mpirun -np 12 pmemd.MPI -O -i heat.mdin -c 6pti.min.rst7 -p 6pti.mod0.parm7 -ref 6pti.min.rst7 -cpin 6pti.cpin -o 6pti.heat.mdout -r 6pti.heat.rst7 -x 6pti.heat.nc

Equilibrating The System

Explicit solvent constant pH molecular dynamics &cntrl imin=0, irest=1, ntx=5, ntpr=1000, ntwx=1000, nstlim=1000000, dt=0.002, ntt=3, tempi=300, temp0=300, gamma_ln=5.0, ig=-1 ntp=1, ntc=2, ntf=2, cut=8, ntb=2, iwrap=1,ioutfm=1, /

mpirun -np 12 pmemd.MPI -O -i equi.mdin -p 6pti.mod0.parm7 -c 6pti.heat.rst7 -cpin 6pti.cpin -o 6pti.equil.mdout -cpout 6pti.equil.cpout -r 6pti.equil.rst7 -x 6pti.equil.nc -cprestrt 6pti.equil.cpin

heating is NVT, and equilibrating is NPT, "Another equilibration issue involves the use of the NPT ensemble. If the system is not reasonably close to a good minimum when you initiate a constant pressure simulation, you can get P-V oscillations. One way to get a decent starting structure for NPT is to start with a short NVT simulation. This is what we do in the Amber water demo".   --George Seibel

It is better use CPU code than GPU code for equilibration,since GPU code does not automatically reorganize grid cells and if the density changed too much, GPU code will fail.

Production Run

Explicit solvent constant pH molecular dynamics &cntrl imin=0, irest=1, ntx=5,ntxo=2, ntpr=1000, ntwx=1000, nstlim=1000000, dt=0.002, ntt=3, tempi=300, temp0=300, gamma_ln=5.0, ig=-1, ntc=2, ntf=2, cut=8, iwrap=1, ioutfm=1, icnstph=2, ntcnstph=100, solvph=3.0, ntrelax=100, saltcon=0.1, /

*icnstph=1 implicit
=2 explicit

for i in $(seq 0 14) ; do sed "s/solvph=3.0/solvph=$i/" md.mdin > pH_$i.mdin; done

for ph in $(seq 0 14); do mpirun -np 2 pmemd.MPI -O -i pH_${ph}.mdin -p 6pti.mod0.parm7 -c 6pti.equil.rst7 -cpin 6pti.cpin -o ${ph}/6pti.md0.mdout -cpout ${ph}/6pti.md0.cpout -r ${ph}/6pti.md0.rst7 -x ${ph}/6pti.md0.nc -cprestrt ${ph}/6pti.md0.cpin & done

Calculate pKas

for ph in $(seq 0 14); do cphstats -i 6pti.cpin ./${ph}/6pti.md0.cpout -o pH${ph}_calcpka.dat --population pH${ph}_populations.dat; done

This will generate 2 files for each pH: pH<ph>_calcpka.dat & pH<ph>_populations.dat

amber_pka.pl 0_14

Multiple Run

 

 

 

 

 

 

 

 

 

Reference:

http://jswails.wikidot.com/explicit-solvent-constant-ph-md

http://dev-archive.ambermd.org/201402/0118.html

http://archive.ambermd.org/201508/0504.html

http://ambermd.org/Questions/equilibration.html