Research:amber modification
nint in sande needs to be replaced by dinit in pmemmd
Modifying Sander (amber 8 by Darden, in house) to accept coupling parameter (Chona)
Note: The instructions here refer to the source directory in lela.tacc.utexas.edu. Except for changing <AMBER_TOP>, these instructions should be platform-independent.
- Go to <AMBER_TOP>/SRC/sander. The sander source files are in ~chona/amoeba/amber06a7l5_p/SRC/sander, so <AMBER_TOP>=~chona/amoeba/amber06a7l5_p.
- Edit file called amoeba_mdin.f. This file contains the source to the amoeba_mdin module, which contains the variables from the input file.
- Add the parameter(s) to the variable definitions:
lela% more amoeba_mdin.f
...
module amoeba_mdin
implicit none
private
integer,save :: do_amoeba=0,do_amoeba_valence=1,do_amoeba_nonbond=1, &
do_vdw_taper=1,do_vdw_longrange=1
_REAL_,save :: sor_coefficient = 0.75d0
_REAL_,save :: dipole_scf_tol = 0.01d0
integer,save :: dipole_scf_iter_max = 50
_REAL_,save :: ee_dsum_cut=7.d0
_REAL_,save :: ee_damped_cut=4.5d0
_REAL_,save :: vdw_taper = 0.9d0
_REAL_,save :: thole_expon_coeff=0.39d0
_REAL_,save :: compress = 0.000046d0
_REAL_,save :: lamda = 0.0
logical, save :: verbose=.false.
integer,save :: beeman_integrator = 0
public AMOEBA_read_mdin,do_amoeba,do_amoeba_valence, &
do_amoeba_nonbond,verbose,beeman_integrator, &
sor_coefficient,dipole_scf_tol,dipole_scf_iter_max, &
ee_dsum_cut,ee_damped_cut,thole_expon_coeff,vdw_taper, &
compress,do_vdw_taper,do_vdw_longrange,lamda
- Include the new parameter in the amoeba_mdin namelist
lela% more amoeba_mdin.f
namelist/amoeba/do_amoeba,do_amoeba_valence,do_amoeba_nonbond, &
do_bond,do_ureyb,do_reg_angle, &
do_trig_angle,do_opbend,do_torsion,do_str_torsion, &
do_pi_torsion,do_strbend,do_torsion_torsion, &
do_recip,do_adjust,do_direct,do_self, &
do_vdw,do_induced,verbose,beeman_integrator, &
sor_coefficient,dipole_scf_tol, &
dipole_scf_iter_max,ee_dsum_cut,ee_damped_cut, &
thole_expon_coeff,vdw_taper,do_vdw_taper,do_vdw_longrange, &
compress, &
lamda
- Include the variable in amoeba_vdw.f
lela% more amoeba_mdin.f
subroutine AM_VDW_DIRECT_ene_frc_i(i,ipairs,numtot,xk,yk,zk, &
crd,ene_vdw,frc,virial)
use nblist, only: bckptr,imagcrds,tranvec
use constants, only : zero,one,two,three,four,five,seven
use amoeba_mdin, only : do_vdw_taper, lamda
integer,intent(in) :: i,ipairs(*),numtot
_REAL_,intent(in) :: xk,yk,zk,crd(3,*)
_REAL_,intent(inout) :: ene_vdw,frc(3,*),virial(3,3)
- Include the new parameter in the input file
lela% more vdw1.mdin &amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba=1,do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=1, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1, do_vdw=1,verbose=.false.,do_vdw_longrange=1,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.05,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, lamda=0.6
Modifying PMEMD (Oscar)
This is an instruction of adding softcore parameters in pmemd.
The original files are in /home/oscar/pmemd/src-org, while the modified files are in /home/oscar/pmemd/src
Four files have been modified: nmr_calls.fpp, mdin_amoeba_dat.fpp, amoeba_direct.fpp, amoeba_induced.fpp.
You can download these new files with softcore modification from http://biomol.bme.utexas.edu/~oscar/download/pmemd/src/
Just copy the four files, overwrite your old ones and compile.
Old files can be found in http://biomol.bme.utexas.edu/~oscar/download/pmemd/src-org/ Just in case you want to change back to original version.
- In order to use restraint, some change should be made in nmr_calls.fpp file by deleting some lines in subrountine restal
#ifdef AMOEBA
! For Amoeba we don't currently support file redirections because we don't
! understand the implications.
if (iamoeba .ne. 0) then
do m = 1, 8
if (iredir(m) .ne. 0) then
write(mdout, *) 'No Amoeba support for file redirections!'
call mexit(6, 1)
end if
end do
! Check to be sure we did not read a card for an nmr option that has been
! dropped in pmemd. These include NOESY, SHIFTS, PCSHIFT and DIPOLE.
! Note that the indices checked are hardwired - what a kludge!
else if (iredir(4) .ne. 0 .or. iredir(5) .ne. 0 .or. iredir(7) .ne. 0 .or. &
iredir(8) .ne. 0) then
write(mdout,*) &
'No support for NOESY, SHIFTS, PCSHIFT or DIPOLE file redirections!'
write(mdout,*)'Please use SANDER instead.'
call mexit(6, 1)
end if
#else
! Check to be sure we did not read a card for an nmr option that has been
! dropped in pmemd. These include NOESY, SHIFTS, PCSHIFT and DIPOLE.
! Note that the indices checked are hardwired - what a kludge!
if (iredir(4) .ne. 0 .or. iredir(5) .ne. 0 .or. iredir(7) .ne. 0 .or. &
iredir(8) .ne. 0) then
write(mdout,*) &
'No PMEMD support for NOESY, SHIFTS, PCSHIFT or DIPOLE options!'
write(mdout,*)'Please use SANDER instead.'
call mexit(6, 1)
end if
#endif /* AMOEBA */
And
#ifdef AMOEBA
! Amoeba does not currently support stuff found in &rst...
if (iamoeba .ne. 0) then
write(mdout, *) 'No Amoeba support for &rst namelist-based options!'
call mexit(6, 1)
end if
#endif /* AMOEBA */
Basically everything between "#indef AMOEBA" and "#endif /*AMOEBA*/" should be deleted.
- Add new variables in the mdin name list. Edit the file mdin_amoeba_dat.fpp.
First, increase the number of variables, three new integers and two new doubles.
Second, add new variables to the end of namelists including both interger and double types
mdin_amoeba_int_cnt = 26
integer do_amoeba_valence, do_amoeba_nonbond, do_bond, &
do_ureyb, do_reg_angle, do_trig_angle, do_opbend, &
do_torsion, do_pi_torsion, do_strbend, &
do_torsion_torsion, do_str_torsion, do_recip, &
do_adjust, do_direct, do_self, do_vdw, &
do_induced, do_vdw_taper, do_vdw_longrange, &
beeman_integrator, dipole_scf_iter_max, &
amoeba_verbose, soft_expo, lig_idx1, lig_idx2
common / mdin_amoeba_int / do_amoeba_valence, do_amoeba_nonbond, do_bond, &
do_ureyb, do_reg_angle, do_trig_angle, do_opbend, &
do_torsion, do_pi_torsion, do_strbend, &
do_torsion_torsion, do_str_torsion, do_recip, &
do_adjust, do_direct, do_self, do_vdw, &
do_induced, do_vdw_taper, do_vdw_longrange, &
beeman_integrator, dipole_scf_iter_max, &
amoeba_verbose,&
soft_expo, lig_idx1, lig_idx2
integer, parameter :: mdin_amoeba_dbl_cnt = 9
double precision compress, dipole_scf_tol, ee_dsum_cut, &
ee_damped_cut, sor_coefficient, &
thole_expon_coeff, vdw_taper, &
soft_lamda, soft_alpha
common / mdin_amoeba_dbl / compress, dipole_scf_tol, ee_dsum_cut, &
ee_damped_cut, sor_coefficient, &
thole_expon_coeff, vdw_taper, &
soft_lamda, soft_alpha
soft_expo: exponential of softcore equation (old 4; new 5)
soft_alpha: coefficient of softcore (old 0.5, new 0.7)
soft_lamda: scaling factor of decoupling usually changing from 1.0 to 0.0 linearly
lig_idx1,lig_idx2: first and last index of the ligand in the pdb file of the protein holo structure
- Add those variables in the program that related to softcore calculation, amoeba_direct.fpp
In amoeba_direct.fpp, add those variables in subroutine am_vdw_direct_ene_frc_i
use mdin_amoeba_dat_mod, only : do_vdw_taper, soft_lamda, soft_expo, soft_alpha, lig_idx1, lig_idx2
- AMBER exclude 1-2, 1-3 intramolecular interactions and pairs within the same polarization group here so no effect on intramolecular interactions
- Inplement softcore calculation in amoeba_direct.fpp
In the subroutine am_vdw_direct_ene_frc_i, add the following code
if((atm_i .ge. lig_idx1 .and. atm_i .le. lig_idx2) .or. (atm_j .ge. lig_idx1 .and. atm_j .le. lig_idx2)) then eps = eps * soft_lamda ** soft_expo t1 = (1.d0 + vdw_buf_delta)**7 / (soft_alpha * (1.d0 - soft_lamda)**2 + (rho + vdw_buf_delta)**7) t2 = (1.d0 + vdw_buf_gamma) / (soft_alpha * (1.d0 - soft_lamda)**2 + rho**7 + vdw_buf_gamma) dt1drho = (-7.d0 * (rho + vdw_buf_delta)**6 * t1) / (soft_alpha * (1.d0 - soft_lamda)**2 + (rho + vdw_buf_delta)**7) dt2drho = (-7.d0 * rho**6 * t2) / (soft_alpha * (1.d0 - soft_lamda)**2 + rho**7 + vdw_buf_gamma) endif
- Change the damping factor of metal ions, such as Ca2+ and Mg2+
In amoeba_induced.fpp, insert the following lines in subroutine amoeba_induced_dat_final_setup
if ( abs(atm_mass(n) - 40.078) < 0.0001d0 ) & write (6,*) "Ca damp",n,atm_mass(n),polarizability(n) if ( abs(atm_mass(n) - 40.078) < 0.0001d0 ) & atm_qterm(n) = 1 /1.35d0**3*atm_qterm(n) if ( abs(atm_mass(n) - 24.305) < 0.0001d0 ) & write (6,*) "Mg damp",n,atm_mass(n),polarizability(n) if ( abs(atm_mass(n) - 24.305) < 0.0001d0 ) & atm_qterm(n) = 1 /1.60d0**3*atm_qterm(n)
- Do make clean and then make to update the modifacatiions
mdin for amber 8
short md, nve ensemble &cntrl ntx=7, irest=1, nstlim=300000, ntpr=500, ntwr=500,ntave=500, nscm=500,ntwx=500, dt=0.001, vlimit=10.0, cut=9.,maxcyc=50,ntmin=2,imin=0, ntt=1, temp0=298.0,tempi=0.0,tautp=0.5, nmropt=1, / &ewald nfft1=64,nfft2=64,nfft3=64, skinnb=0.8,nbtell=0,order=5,ew_coeff=0.45, / &amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba=1,do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=1, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1, do_vdw=1,verbose=.false.,do_vdw_longrange=1,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.05,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, / &wt type='END' &end LISTOUT=POUT DISANG=bnd.rst
mdin for PMEMD/amber9
short md, nve ensemble &cntrl ntx=7, irest=1, nstlim=3000000, ntpr=500, ntwr=500,ntave=500, nscm=500,ntwx=500, dt=0.001, vlimit=10.0, cut=9.,maxcyc=50,ntmin=2,imin=0, ntt=1, temp0=298.0,tempi=0.0,tautp=0.5, nmropt=1, iamoeba=1, / &ewald nfft1=64,nfft2=64,nfft3=64, nbtell=0,order=5,ew_coeff=0.45, / &amoeba do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1,do_amoeba_nonbond=1, dipole_scf_tol = 0.01,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, beeman_integrator=1,<b>lig_idx1=3224,lig_idx2=3241,</b> <b>soft_alpha=0.7,soft_expo=5,soft_lamda=0.6</b> / &wt type='END' &end LISTOUT=POUT DISANG=bnd.rst
The last two lines are for distance restraint.
nfft should be 20% bigger than box length. Usable values:
data multi / 2, 3, 4, 5, 6, 8, 9, 10, 12, 15,
& 16, 18, 20, 24, 25, 27, 30, 32, 36, 40,
& 45, 48, 50, 54, 60, 64, 72, 75, 80, 81,
& 90, 96, 100, 108, 120, 125, 128, 135, 144, 150,
& 160, 162, 180, 192, 200, 216, 225, 240, 243, 250,
& 256, 270, 288, 300, 320, 324, 360, 375, 384, 400,
& 405, 432, 450, 480, 486, 500, 512, 540, 576, 600 /
For npt, add this line (see amber manual. ntb=2 is pressure control. ntp=1 isoropic, pres0 is target) ntb=2, ntp=1, pres0=1.0, taup=50.0 For position restraint: ntr=1,restraint_wt=10.0,restraintmask=':1-224'
To run pmemd:
creat a host file, containing the following lines node25:4 node24:4 run mpdboot, not the ncpus should be less than the total avail in host file. /opt/mpich2-fc7/bin/mpdboot -n 1 -f ./host --ncpus=4 run pmemed: /opt/mpich2-fc7/bin/mpirun -n 4 ~/pmemd/src/pmemd -i mdin_pmemd -c ele1.inpcrd -p ele1.prmtop -o mdout.pmemd -r restrt.pmemd -x mdcrd.pmemd 2> /dev/full &
Modifying pre-AMBER9 (Oscar)
This is an instruction of adding softcore parameters in AMBER 8.
The original files are in /home/oscar/DARDEN_amber/SRC/sander/, while the modified files are in /home/oscar/data2-oscar/amber_mod/SRC/sander/
/home/oscar/data2-oscar/amber-atm-soft: two atoms specified for soft-core
/home/oscar/data2-oscar/amber-range-soft: a range of (ligand atoms) is specified for soft-core
Two files have been modified: amoeba_mdin.f and amoeba_vdw.f.
You can also download these new files with softcore modification from http://biomol.bme.utexas.edu/~oscar/download/amber8/src-mod/
Just copy the two files, overwrite your old ones and compile.
Old files can still be found in http://biomol.bme.utexas.edu/~oscar/download/amber8/src-org/ Just in case you want to change back to original version.
- Add the parameter(s) to the variable definitions in amoeba_mdin.f and assign initial values to them:
...
module amoeba_mdin
implicit none
private
integer,save :: do_amoeba=0,do_amoeba_valence=1,do_amoeba_nonbond=1, &
do_vdw_taper=1,do_vdw_longrange=1
_REAL_,save :: sor_coefficient = 0.75d0
_REAL_,save :: dipole_scf_tol = 0.01d0
integer,save :: dipole_scf_iter_max = 50
_REAL_,save :: ee_dsum_cut=7.d0
_REAL_,save :: ee_damped_cut=4.5d0
_REAL_,save :: vdw_taper = 0.9d0
_REAL_,save :: thole_expon_coeff=0.39d0
_REAL_,save :: compress = 0.000046d0
_REAL_,save :: softcore_lamda = 1.0d0
_REAL_,save :: softcore_alpha = 0.5d0
_REAL_,save :: softcore_expo = 4
_REAL_,save :: lig_idx1 = 1
_REAL_,save :: lig_idx2 = 1
logical, save :: verbose=.false.
integer,save :: beeman_integrator = 0
public AMOEBA_read_mdin,do_amoeba,do_amoeba_valence, &
do_amoeba_nonbond,verbose,beeman_integrator, &
sor_coefficient,dipole_scf_tol,dipole_scf_iter_max, &
ee_dsum_cut,ee_damped_cut,thole_expon_coeff,vdw_taper, &
compress,do_vdw_taper,do_vdw_longrange, &
softcore_lamda,softcore_alpha,softcore_expo, &
lig_idx1, lig_idx2
softcore_expo: exponential of softcore equation with initial value 4
softcore_alpha: coefficient of softcore with initial value 0.5d0
softcore_lamda: scaling factor of decoupling usually changing from 1.0 to 0.0 linearly with initial value 1.0d0
lig_idx1,lig_idx2: first and last index of the ligand in the pdb file of the protein holo structure with initial values 1
- Include the new parameters in the amoeba_mdin namelist
namelist/amoeba/do_amoeba,do_amoeba_valence,do_amoeba_nonbond, &
do_bond,do_ureyb,do_reg_angle, &
do_trig_angle,do_opbend,do_torsion,do_str_torsion, &
do_pi_torsion,do_strbend,do_torsion_torsion, &
do_recip,do_adjust,do_direct,do_self, &
do_vdw,do_induced,verbose,beeman_integrator, &
sor_coefficient,dipole_scf_tol, &
dipole_scf_iter_max,ee_dsum_cut,ee_damped_cut, &
thole_expon_coeff,vdw_taper,do_vdw_taper,do_vdw_longrange, &
compress, softcore_lamda,softcore_alpha,softcore_expo, &
lig_idx1,lig_idx2
- Include the variables in amoeba_vdw.f
subroutine AM_VDW_DIRECT_ene_frc_i(i,ipairs,numtot,xk,yk,zk, &
crd,ene_vdw,frc,virial)
use nblist, only: bckptr,imagcrds,tranvec
use constants, only : zero,one,two,three,four,five,seven
use amoeba_mdin, only : do_vdw_taper, lig_idx1,lig_idx2, &
softcore_lamda,softcore_alpha,softcore_expo
integer,intent(in) :: i,ipairs(*),numtot
_REAL_,intent(in) :: xk,yk,zk,crd(3,*)
_REAL_,intent(inout) :: ene_vdw,frc(3,*),virial(3,3)
- Insert the softcore algorithm with new parameters in the same subroutine
dt1drho = -seven*t1 / (rho + vdw_buf_delta)
dt2drho = -seven*t2 * (rho6 / (rho7 + vdw_buf_gamma))
drhodr = one / rad
if ((lig_idx1 <= i .and. i <= lig_idx2) .or. (lig_idx1 <= j .and. j <= lig_idx2) )then
eps = eps * softcore_lamda ** softcore_expo
t1 = (one + vdw_buf_delta)**7 / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7)
t2 = (one + vdw_buf_gamma) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma)
dt1drho = (-seven * (rho + vdw_buf_delta)**6 * t1) / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7)
dt2drho = (-seven * rho6 * t2) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma)
endif
f = eps*t1*(t2 - two)
dfdr = eps*(dt1drho*(t2 - two) + t1*dt2drho)*drhodr
- Include the new parameter in the input file (mdin)
&amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba=1,do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=1, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1, do_vdw=1,verbose=.false.,do_vdw_longrange=1,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.05,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, softcore_lamda=0.9,softcore_alpha=0.7,softcore_expo=5, lig_idx1=3224,lig_idx2=3241
Here we take lambda=0.9 alpha=0.7 exponent=5 as example. The ligand atom index ranges from 3224 to 3241 (benzamidine in trypsin).
Modifying Sander/AMBER9 (Yue)
Adding softcore parameters in Sander
The modified files are in
/home/yues/amber9/src-range-soft/ (a range of ligand atoms is specified for soft-core) /home/yues/amber9/src-atm-soft/ (two atoms specified for soft-core)
Two files have been modified: amoeba_mdin.f and amoeba_vdw.f.
- Add the parameter(s) to the variable definitions in amoeba_mdin.f and assign initial values to them:
...
module amoeba_mdin
implicit none
private
integer,save :: iamoeba=0,do_amoeba_valence=1,do_amoeba_nonbond=1, &
do_vdw_taper=1,do_vdw_longrange=1
_REAL_,save :: sor_coefficient = 0.75d0
_REAL_,save :: dipole_scf_tol = 0.01d0
integer,save :: dipole_scf_iter_max = 50
_REAL_,save :: ee_dsum_cut=7.d0
_REAL_,save :: ee_damped_cut=4.5d0
_REAL_,save :: vdw_taper = 0.9d0
_REAL_,save :: thole_expon_coeff=0.39d0
_REAL_,save :: compress = 0.000046d0
_REAL_,save :: softcore_lamda = 1.0d0
_REAL_,save :: softcore_alpha = 0.5d0
_REAL_,save :: softcore_expo = 4
'''_REAL_,save :: lig_idx1 = 1'''
_REAL_,save :: lig_idx2 = 1</span><br> integer, save :: amoeba_verbose = 0<br> integer,save :: beeman_integrator = 0<br> public AMOEBA_read_mdin,iamoeba,do_amoeba_valence, &<br> do_amoeba_nonbond,amoeba_verbose,beeman_integrator, &<br> sor_coefficient,dipole_scf_tol,dipole_scf_iter_max, &<br> ee_dsum_cut,ee_damped_cut,thole_expon_coeff,vdw_taper, &<br> compress,<span style="color: rgb(255, 0, 0);">do_vdw_taper,do_vdw_longrange, &
softcore_lamda,softcore_alpha,softcore_expo, &
lig_idx1, lig_idx2</span><br><br>
softcore_expo: exponential of softcore equation with initial value 4
softcore_alpha: coefficient of softcore with initial value 0.5d0
softcore_lamda: scaling factor of decoupling usually changing from 1.0 to 0.0 linearly with initial value 1.0d0
lig_idx1,lig_idx2: first and last index of the ligand in the pdb file of the protein holo structure with initial values 1
- Include the noew parameters in the subroutine AMOEBA_read_mdin(nf) namelist
namelist/amoeba/do_amoeba_valence,do_amoeba_nonbond, & do_bond,do_ureyb,do_reg_angle, & do_trig_angle,do_opbend,do_torsion,do_str_torsion, & do_pi_torsion,do_strbend,do_torsion_torsion, & do_recip,do_adjust,do_direct,do_self, & do_vdw,do_induced,amoeba_verbose,beeman_integrator, & sor_coefficient,dipole_scf_tol, & dipole_scf_iter_max,ee_dsum_cut,ee_damped_cut, & thole_expon_coeff,vdw_taper,do_vdw_taper,do_vdw_longrange, & compress, '''''softcore_lamda,softcore_alpha,softcore_expo, & lig_idx1, lig_idx2'''''
- Include the variables in amoeba_vdw.f
subroutine AM_VDW_DIRECT_ene_frc_i(i,ipairs,numtot,xk,yk,zk, & crd,ene_vdw,frc,virial) use nblist, only: bckptr,imagcrds,tranvec use constants, only : zero,one,two,three,four,five,seven use amoeba_mdin, only : do_vdw_taper,'''''lig_idx1,lig_idx2, & softcore_lamda,softcore_alpha,softcore_expo''''' integer,intent(in) :: i,ipairs(*),numtot _REAL_,intent(in) :: xk,yk,zk,crd(3,*) _REAL_,intent(inout) :: ene_vdw,frc(3,*),virial(3,3)
- Insert the softcore algorithm with new parameters in the same subroutine
/home/yues/amber9/src-range-soft/ (a range of ligand atoms is specified for soft-core)
if ( delr2 < vdw_switch_off_2 )then jt = vdw_atom_type(j) eps = vdw_epsilon(jt,it) rad = vdw_radius(jt,it) delr = sqrt(delr2) rho = delr / rad rho6 = rho**6 rho7 = rho6*rho t1 = ((one + vdw_buf_delta) / (rho + vdw_buf_delta))**7 t2 = (one + vdw_buf_gamma) / (rho7 + vdw_buf_gamma) dt1drho = -seven*t1 / (rho + vdw_buf_delta) dt2drho = -seven*t2 * (rho6 / (rho7 + vdw_buf_gamma)) drhodr = one / rad ''''' if ((lig_idx1 <= i .and. i <= lig_idx2) .or. (lig_idx1 <= j .and. j <= lig_idx2) )then eps = eps * softcore_lamda ** softcore_expo t1 = (one + vdw_buf_delta)**7 / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7) t2 = (one + vdw_buf_gamma) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma) dt1drho = (-seven * (rho + vdw_buf_delta)**6 * t1) / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7) dt2drho = (-seven * rho6 * t2) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma) endif ''''' f = eps*t1*(t2 - two) dfdr = eps*(dt1drho*(t2 - two) + t1*dt2drho)*drhodr
/home/yues/amber9/src-atm-soft/ (two atoms specified for soft-core)
if ( delr2 < vdw_switch_off_2 )then jt = vdw_atom_type(j) eps = vdw_epsilon(jt,it) rad = vdw_radius(jt,it) delr = sqrt(delr2) rho = delr / rad rho6 = rho**6 rho7 = rho6*rho t1 = ((one + vdw_buf_delta) / (rho + vdw_buf_delta))**7 t2 = (one + vdw_buf_gamma) / (rho7 + vdw_buf_gamma) dt1drho = -seven*t1 / (rho + vdw_buf_delta) dt2drho = -seven*t2 * (rho6 / (rho7 + vdw_buf_gamma)) drhodr = one / rad ''''' if ((lig_idx1 == i .and. i == lig_idx2) .or. (lig_idx1 == j .and. j == lig_idx2) )then eps = eps * softcore_lamda ** softcore_expo t1 = (one + vdw_buf_delta)**7 / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7) t2 = (one + vdw_buf_gamma) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma) dt1drho = (-seven * (rho + vdw_buf_delta)**6 * t1) / (softcore_alpha * (1 - softcore_lamda)**2 + (rho + vdw_buf_delta)**7) dt2drho = (-seven * rho6 * t2) / (softcore_alpha * (1 - softcore_lamda)**2 + rho7 + vdw_buf_gamma) endif '''''f = eps*t1*(t2 - two) dfdr = eps*(dt1drho*(t2 - two) + t1*dt2drho)*drhodr
- Include the new parameter in the input file (mdin)
&amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=1, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1, do_vdw=1,amoeba_verbose=0,do_vdw_longrange=1,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.05,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, '''softcore_lamda=0.4,softcore_alpha=0.7,softcore_expo=5, lig_idx1=12,lig_idx2=13 '''
Changing the damping factor of metal ions,such as Ca2+ and Mg2+
There are three files modified: amoeba_interface.f, amoeba_induced.f and sander.f.
- In amoeba_induced.f, add mass variable in function AM_INDUCED_readparm.
function AM_INDUCED_readparm(nf,num_atoms,'''''mass''''') integer :: AM_INDUCED_readparm integer,intent(in) :: nf,num_atoms '''''_REAL_ mass(*)'''''
- Then insert the following lines in the same function
do n = 1,num_atoms if ( polarizability(n) > polmin )then is_polarizable(n) = .true. sq_polinv(n) = 1.d0 / sqrt(polarizability(n)) '''''if ( abs(mass(n) - 40.078) < 0.0001d0 ) & write (6,*) "Ca damp",n,mass(n),polarizability(n) if ( abs(mass(n) - 40.078) < 0.0001d0 ) & sq_polinv(n) = 1 /1.35d0**3*sq_polinv(n) if ( abs(mass(n) - 24.305) < 0.0001d0 ) & write (6,*) "Mg damp",n,mass(n),polarizability(n) if ( abs(mass(n) - 24.305) < 0.0001d0 ) & sq_polinv(n) = 1 /1.60d0**3*sq_polinv(n)''''' endif enddo AM_INDUCED_readparm = 1 do_flag = ibset(do_flag,VALID_BIT) end function AM_INDUCED_readparm
- After that, add mass variables and its definition in any programs or subroutines that is related to this function.
In amoeba_interface.f
subroutine AM_NONBOND_readparm(nf,natom,'''''mass''''') use amoeba_multipoles, only : AM_MPOLE_readparm use amoeba_adjust, only : AM_ADJUST_readparm use amoeba_vdw, only : AM_VDW_read_parm use amoeba_induced, only : AM_INDUCED_readparm use amoeba_recip, only : AM_RECIP_allocate use amoeba_runmd, only : AM_RUNMD_init integer,intent(in) :: nf,natom ''''' _REAL_ mass(*)''''' integer :: mpole_valid,adjust_valid,vdw_valid,polar_valid mpole_valid = AM_MPOLE_readparm(nf,natom) adjust_valid = AM_ADJUST_readparm(nf) vdw_valid = AM_VDW_read_parm(nf) polar_valid = AM_INDUCED_readparm(nf,natom,'''''mass''''')
subroutine AMOEBA_readparm(nf,ntf,ntc,numatoms,'''''mass''''')use amoeba_mdin, only : iamoeba integer,intent(in) :: nf,numatoms integer,intent(inout) :: ntf,ntc !set these to 8,1 if amoeba-friendly prmtop !so regular amber valence & shake not run '''''_REAL_ mass(*) ''''' call AMOEBA_check_parm_legal(nf)if ( iamoeba /= 1 )return
call AM_VAL_readparm(nf,ntf,ntc) call AM_NONBOND_readparm(nf,numatoms,'''''mass''''')end subroutine AMOEBA_readparm
In sander.f
call rdparm2(x,ix,ih,ipairs,8) call AMOEBA_readparm(8,ntf,ntc,natom,'''''x(lmass)''''')! ntf,ntc get reset if amoeba prmtop if (qmmm_nml%ifqnt) then call read_qmmm_namelist_and_allocate(igb, ih, ix, x, cut, use_pme, ntb) end if call mdread2(x,ix,ih,ipairs) ! --- alloc memory for decomp module that needs info from mdread2
PMEMD.AMOEBA Modifications in Amber 12
a. RESPA
- Respa performs multiple time step molecular dynamics using the reversible reference system propagation algorithm (r-RESPA) via a Verlet core with the potential split into fast- and slow-evolving portions.
- Code Modification
- Added a new file (amoeba_force_respa.F90) to calculate the RESPA energy and force
!every nrespa step, calculate the nonbond interactions
if (need_pot_enes) then
call am_nonbond_eval(atm_cnt, crd, frc, img_frc, virial_lcl, &
pot_ene%vdw, pot_ene%elec, pot_ene%polar, &
pot_ene%vdw_14, diprms, dipiter, netfrcs)
end if
...
!Other conditions, calculate the valence interactions only
if (.not. need_pot_enes) then
call am_val_eval(crd, frc, virial_lcl, &
pot_ene%bond, pot_ene%angle, pot_ene%dihedral)
end if
- Important Modifications in runmd.F90
use amoeba_force_respa_mod
…
! When using respa, call bonded force or nonbonded force at different conditions
!bond
need_pot_enes = .false.
#ifdef MPI
call amoeba_force_respa(atm_cnt, crd, frc, gbl_img_atm_map, gbl_atm_img_map, &
my_atm_lst, new_list, need_pot_enes, need_virials, amba_pot_ene,
si(si_diprms), si(si_dipiter), virial)
#else
call amoeba_force_respa(atm_cnt, crd, frc, gbl_img_atm_map, gbl_atm_img_map, &
new_list, need_pot_enes, need_virials, amba_pot_ene, &
si(si_diprms), si(si_dipiter), virial)
#endif /* MPI */
!nonbond
#ifdef MPI
call amoeba_force_respa(atm_cnt, crd, frc, gbl_img_atm_map, gbl_atm_img_ma
my_atm_lst, new_list, need_pot_enes, need_virials, amba_pot_ene,
si(si_diprms), si(si_dipiter), virial)
#else
call amoeba_force_respa(atm_cnt, crd, frc, gbl_img_atm_map, gbl_atm_img_ma
new_list, need_pot_enes, need_virials, amba_pot_ene, &
si(si_diprms), si(si_dipiter), virial)
#endif
!Update nonbond contribution to velocity, every nrespa steps
do j = 1, atm_cnt
wfac = atm_mass_inv(j) * dtx
vel(:, j) = vel(:, j) + frc(:, j) * wfac
if (onstep) then
vel(:,j) = vel(:,j) + frc_nb(:,j) * wfac * nrespa
endif
end do
- Input Example
short md, nve ensemble &cntrl irest=1,ntx=7, nstlim=1000, ntpr=10, ntwr=100,ntave=100, nscm=100,ntwx=100, dt=0.00025, nrespa=10, vlimit=10.0, cut=12.0,maxcyc=50,ntmin=2,imin=0, ntt=4, temp0=298.0,tempi=298.0,tautp=0.5, iamoeba=1 / &ewald nfft1=36,nfft2=36,nfft3=36, skinnb=2.0,nbtell=0,order=5,ew_coeff=0.45, / &amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=0, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1,amoeba_verbose=0, do_vdw=1,do_vdw_longrange=0,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.0001,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=8, /
dt=0.00025 Time step of valence interactions
nrespa=10 Update nonbonded interactions every nrespa steps
b. Bussi
- An extension of the Berendsen thermostat in which a random force is added to ensure the correct kinetic energy distribution
- Modifications in runmd.F90
if (ntt .eq. 4 .and. onstep) then
temp = eke / fac(1)
if (temp .eq. 0.0d0) temp = 0.1d0
c = exp(-dt * nrespa/tautp)
d = (1.0d0-c) * (temp0/temp) / (rndfp + rndfs)
call gauss(0.d0,1.d0,r)
s = 0.0d0
do i = 1, (rndfp + rndfs-1)
call gauss(0.d0,1.d0,si_bussi)
s = s + si_bussi*si_bussi
end do
scaltp = c + (s+r*r)*d + 2.0d0*r*sqrt(c*d)
scaltp = sqrt(scaltp)
if (r+sqrt(c/d) .lt. 0.0d0) scaltp = -scaltp
#ifdef MPI
do atm_lst_idx = 1, my_atm_cnt
j = my_atm_lst(atm_lst_idx)
#else
do j = 1, atm_cnt
#endif
vel(:, j) = vel(:, j) * scaltp
end do
end if !(bussi)
- Input Example
short md, nve ensemble &cntrl irest=1,ntx=7, nstlim=1000, ntpr=10, ntwr=100,ntave=100, nscm=100,ntwx=100, dt=0.00025, vlimit=10.0,nrespa=10, cut=12.0,maxcyc=50,ntmin=2,imin=0, ntt=4, temp0=298.0,tempi=298.0,tautp=0.5, iamoeba=1 / &ewald nfft1=36,nfft2=36,nfft3=36, skinnb=2.0,nbtell=0,order=5,ew_coeff=0.45, / &amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=0, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1,amoeba_verbose=0, do_vdw=1,do_vdw_longrange=0,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.0001,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=8, /
ntt=4 Enable bussi thermostat
c. REMD
- Copy remd_exchg.F90 and remd.F90 over from regular PMEMD
- Add REMD loop to runmd.F90 and amoeba_runmd.F90 (similar to regular PMEMD)
- Other files with minor changes: bintraj.F90, file_io_dat.F90, get_cmdline.F90, master_setup.F90, mdin_ctrl_dat.F90, multipmemd.F90, parallel_dat.F90, pmemd.amoeba.F90, runfiles.F90
- Input Example
short md, nve ensemble &cntrl irest=0, nstlim=4000,numexchg=1000, ntpr=4000,ntwr=4000,ntave=4000, nscm=4000,ntwx=4000, dt=0.00025, vlimit=10.0, nrespa=10 cut=12.0,maxcyc=50,ntmin=2,imin=0, ntt=4, temp0=XXXXX, ig=RANDOM_NUMBER, tautp=1.0, iamoeba=1, / &ewald nfft1=40,nfft2=40,nfft3=40, skinnb=0.8,nbtell=0,order=5,ew_coeff=0.45, / &amoeba do_induced=1,do_recip=1,do_direct=1,do_adjust=1, do_amoeba_nonbond=1,do_amoeba_valence=1,beeman_integrator=0, do_bond=1,do_ureyb=1,do_reg_angle=1,do_trig_angle=1, do_opbend=1,do_torsion=1,do_str_torsion=1,do_pi_torsion=1,do_strbend=1, do_torsion_torsion=1,amoeba_verbose=0, do_vdw=1,do_vdw_longrange=1,do_vdw_taper=1, do_self=1,dipole_scf_tol = 0.01,dipole_scf_iter_max=30, sor_coefficient=0.6,ee_damped_cut=4.5,ee_dsum_cut=6.8, /
nstlim: Number of steps to run in each exchange numexchg: Number of exchanges