Research ex:AMOEBA in Amber
sander and pmemd modification for ligand binding free enery calculations
Important: read me first
Now (as of June 2013) you can use the latest tinker (v6) and parameter files to set up the system in TINKER and then convert for amber run.
- If you have an spherical ion (e.g. Na+, Cl-, Mg++in the system), you may need to modify the analout (see below) to replace multipole type "None" with Z-then-X
- If you have an octahedron box, please modify the prmtop and inpcrd generated by tinker_to_amber manually according to Quick guide below
(will fix in the future tinker_to_amber.F90).
- You can still use tinker v4.3 for now and download tinker 4.3 here http://dasher.wustl.edu/tinker/
- The difference between v4.3 and v6 parameter format is explained here explained here .
- You can use amoeba_conv.pl to convert v5/6 key files to 4.3 compatible. amoeba_conv.pl:http://biomol.bme.utexas.edu/~yues/download/amoeba_conv.pl
- Distance restraint is not working with the current pmemd.amoeba (see /home/pren/pmemd-oscar/readme)
For proteins and many small molecules, parameters in tinker 4 or 6 compatible format can be found here
For a new small molecule, use POLTYPE to obtained tinker structure and parameter file.
Once you have the tinker files ready, the quick guide section below explains how to convert tinker files to amber by using tinker and amber commands.
A guide to run AMOEBA simulation using sander or pmemd
Executables (for internal use)
Oscar's mod sander in /home/oscar/data2-oscar mod pmemd in /home/oscar/pmemd/src
example mdin and run script files in /home/pren/AMOEBA/solvation/compare/
home/yues/program/amber12/bin/pmemd.amoeba
A Quick Guide
To run sander:
If you start with a PDB, pdbxyz will produce the xyz file with the matching parameters of your choice. But this PDB will not match the xyz file and can not be used for amber simulation so you will need do step 1 below.
1. /opt/tinker/bin/xyzpdb ele.xyz
This produce PDB file from tinker xyz but the atoms have to be in the same order as xyz which is not always true. Use this script fix the order issue: http://biomol.bme.utexas.edu/~yues/download/pdb_crct.pl The most recent tinker analyze (v6.0) also has this option and will output in the v4.3 format recognized by the amoeba_parm below.
2. /opt/tinker/bin/analyze box30 PC > ele.analout #done in tinker v4.3
3. /opt/amber/bin/amoeba_parm -name ele -title "water" > ele.amoeba_parm.out
#this is a program by Tom Darden; Linux executable is at http://biomol.bme.utexas.edu/~pren/amoeba_param/amoeba_parm
4. /opt/amber/bin/sander -i ele.mdin -c ele.inpcrd -p ele.prmtop -o mdout ;
To run pmemd:
The first 3 steps same as above.
4. /opt/mpich2-fc7/bin/mpdboot -n 1 -f ./host --ncpus=4
5. /opt/mpich2-fc7/bin/mpirun -n 4 /home/pren/pmemd-oscar/src/pmemd -i ele.mdin -c ele.inpcrd -p ele.prmtop -o mdout.nvt.pmemd -r restrt.pmemd 2>/dev/null &
A couple of issues about octahedron PBC: TINKER AND AMBER use different oct definition. So the density in amber will be very low once you convert. The conversion script does not recognize octahedron keyword from tinker. See this for how to modify the amber files manually: http://biomol.bme.utexas.edu/wiki/index.php/Research:Amoeba_amber_free#Octahedron
example mdin
minimization using sander
Titlexxxx &cntrl ntx=7, irest=1, nstlim=1000000, ntpr=500, ntwr=500,ntave=500, nscm=500,ntwx=500, dt=0.001, vlimit=10.0, cut=9.,maxcyc=1000,ntmin=2,imin=5, (imin=5 is used to re-evaluate the trajectory during post-processing, as in BAR) iamoeba=1 (this is new for amber10; previous amber didn't use this) / &ewald nfft1=36,nfft2=36,nfft3=36, 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=1, (note do_amoeba=1 is obsolete in amber10) 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, (verbose=false is obsolete) do_self=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, /
NVT using pmemd
&cntrl ntx=1,irest=0, nstlim=1000000, ntpr=100, ntwr=100,ntave=100, 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, iamoeba=1 / &ewald nfft1=36,nfft2=36,nfft3=36, (this depends on box size. minimum 20% bigger) 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=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.01,dipole_scf_iter_max=30, ee_damped_cut=4.5,ee_dsum_cut=6.8, /
NPT using sander
&cntrl irest=0, nstlim=100000, ntpr=100, ntwr=100,ntave=100, 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, ntb=2, ntp=1, pres0=1.0, taup=1.0 ...
direcories in /home/pren/AMOEBA/solvation/compare/
cube: nvt, cubic water box bysander8
cube-pme: nvt, pmemd 9
sander-nvt: truncated octahedron (TO) long nvt run, sander 8
pmemd-nvt: TO, nvt, pmemd 9
sander-long-npt: truncated octahedron long npt run sander 8
pmemd-test-npt: short npt, TO, pmemd
pmemd-nvt/
check softcore parameters defaults in pmemd (sander has initilized correctly).
i.e. without softcore-line lines in mdin, lig_idx is initialized to 0 so it is fine.
softcore-print usem pmemd from ~/pmemd-sun
Simulation set up using AMBER and amoeba
http://biomol.bme.utexas.edu/~yues/download/amber_auto
The automation tool runs the following steps,
1. Auto parameterization using Poltype. 2. Setting up pro-lig binding system in Tinker. 3. Convert the system from Tinker to Amber format 4. Start Amber simulations (3000 steps minimization, 100ps NPT equilibration, NVT simulation with respa)
Input files needed:
1. protein pdb file 2. ligand sdf file 3. amoebapro4.prm parameter file
To use the script, the following environment variables (with the appropriate paths) need to be exported (same as poltype).
$ export PYTHONPATH=/opt/openbabel-2.3.0/lib/python2.6/site-packages:$PYTHONPATH $ export GDMADIR=/opt/gdma2.2/bin $ export LD_LIBRARY_PATH=/opt/openbabel-2.3.0/lib:$LD_LIBRARY_PATH
$ source ~/.g09.bashrc
Usage
1. setup the pro-lig binding system
./runmd.pl -p protein.pdb -l ligand.sdf -o output_trajectory -n simulation_time(fs)
2. run mpdboot and start amber jobs
./run_amber
Example
/home/yues/automation/new/example
Possible problems.
1. Guassian optimization failed.
We have to start the optimization with a lower basis set(like HF/3-21G) and then optimize with the default basis set (HF/6-31G*).
2. Poltype assigns to many atoms to the same polarize group, and Tinker won't continue.
We need to manually devide the pol group and restart.
For free energy calculations, automation is on the way. The instruction for manual operation is here http://biomol.bme.utexas.edu/wiki/index.php/Research:Amoeba_amber_free
RESPA and Bussi thermostat
Example for amber 10 or 12 (PMEMD only, not in sander) using Bussi and RESPA:
/home/yues/testamber/ele1.0/respa/pmemd12-verlet-respa/
Note ntt=4 triggers Bussi thermostat
short md, nvt ensemble using Bussi thermo and RESPA integrator &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, /
and the REMD example:
/home/yues/program/amber12/test/rem_respa_amoeba/
In AMBER 12, pmemd for amoeba can now use 2 fs time step (without shake or rattle).
The dt is the smallest time step used for integration over valence terms. The nrespa x dt is the time step for nonbond terms. I recommend nrespa=8 for stability.
For RESPA, note "dt=0.00025,nrespa=10", and "beeman_integrator=0" are needed for RESPA. Also the printing parameters, ntpr, ntwr, ntave, nscm, ntwx... have to be an integer multiply of nrespa, otherwise the output will be wrong.
md, nvt, RESPA, Langevin &cntrl
irest=0, nstlim=100, ntpr=200 ntwr=200,ntave=200, nscm=200,ntwx=200, dt=0.00025,nrespa=10, vlimit=10.0, cut=9.0,maxcyc=50,ntmin=2,imin=0, ntt=3, temp0=298.0,tempi=298.0,tautp=0.5, iamoeba=1
/ &ewald
nfft1=36,nfft2=36,nfft3=36, 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=300, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=7,
/
AMBER Softcore modification (w/ LRC)
Information below applied to amber 10 (in progress) Latest amber 9
/home/yue/program/amber9/pmemd-oscar/src_LRC/pmemd
Instruction
soft_atm.txt
see /home/pren/pmemd-oscar/node21-amber12-2013 for an example.
The range of atoms that require softcore vdw modification is stored in a file named "soft_atm.txt". It has to be in the same directory as all the files.
The format of soft_atm.txt:
- These atom numbers should be consistent with xyz file (or pdb file)
- Each line contains only one segment
- If it is a range, use hyphen to connect the first atom and the last atom of the segment
Example of "soft_atm.txt"
45-340 580 789-1007
Note that maximum 20 lines are allowed in the soft_atm.txt file </span>
Softcore keywords
Some keywords related to the vdw softcore should be included in the mdin file(input file of amber)
In both sander10 and pmemd.amba10
soft_lambda=0.9,soft_alpha=0.7,soft_expo=5,vdw_longrange_lambda=0.0
soft_alpha: coefficient of softcore default 0.5d0
soft_lambda: scaling factor of decoupling usually changing from 1.0 to 0.0 default 1.0
soft_expo: exponent of scaling factor default 4
vdw_longrange_lambda: scaling factor of vdw longrange correction, default 1.0 (no LRC); 0.0 when there is full LRC
If you are not doing softcore, no need to add these keywords. Even though you have them, as long as softcore_lamda/soft_lamda is 1, no softcore is done.
If you have non-one value for these softcore lambdas, you must have soft_atm.txt file. Otherwise the job will die and complain the file is missing.
Download
Sander10
The original files are in /home/yues/program/amber10/src_org, while the modified files are in /home/yues/program/amber10/src_LRC or "/home/yues/program/amber10/src/sander"
The patch file (sandermod.patch) and config_amber.h can be found in http://biomol.bme.utexas.edu/~yues/download/sander10.
sander 10 may be much slower than sander 9 if the "!<compile=optimized> " line is missing on the top of amoeba*.f files! This line exist in sander 9 source codes.
Pmemd.amba10
The original files are in /home/yues/program/amber10/src_org/pmemd.amba/, while the modified files are in /home/yues/program/amber10/src/pmemd.amba/src-remd-respa pmemd.amba
OLD: /home/yues/program/amber10/src_LRC/pmemd.amba
The patch file (pmemdmod.patch) and config.h can be found in http://biomol.bme.utexas.edu/~yues/download/pmemd10.
Here is how to patch the source files
- Just put the directory src-org and pmemd.patch at the same directory. You need config.h too.
- Go to src-org, do "make clean",
- Go one level up and do "patch -p0 -i pmemd.patch", src-org will be updated.
- Compile and you are good to go!
Examples and tests
Example files are in: /home/yues/testamber/lrc , for softcore tests
In /home/yues/testamber/zinc, there are tests for zinc damping factor
Sander10 tests and Pmemd.amba10 tests are in sander-10/ and pmemd-10/ respectively.
Modification details
Modifying Sander10
amoeba_mdin.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,am_nbead=1
_REAL_,save :: sor_coefficient = 0.75d0
_REAL_,save :: dipole_scf_tol = 0.01d0
integer,save :: dipole_scf_iter_max = 50
#ifndef DISABLE_AMOEBA_CG
integer,save :: dipole_scf_use_cg = 0
integer,save :: dipole_scf_cg_niter = 4
#endif /* DISABLE_AMOEBA_CG */
_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 :: soft_lambda = 1.0d0
_REAL_,save :: soft_alpha = 0.5d0
_REAL_,save :: soft_expo = 4
_REAL_,save :: vdw_longrange_lambda = 1.0d0
integer, save :: amoeba_verbose = 0
integer,save :: beeman_integrator = 0
!variables and arrays for softcore file
integer,save :: soft_atom_range1(20),soft_atom_range2(20),soft_line
public AMOEBA_read_mdin,iamoeba,do_amoeba_valence, &
do_amoeba_nonbond,amoeba_verbose,beeman_integrator, &
sor_coefficient,dipole_scf_tol,dipole_scf_iter_max, &
#ifndef DISABLE_AMOEBA_CG
dipole_scf_use_cg, dipole_scf_cg_niter, &
#endif /* DISABLE_AMOEBA_CG */
ee_dsum_cut,ee_damped_cut,thole_expon_coeff,vdw_taper, &
compress,do_vdw_taper,do_vdw_longrange,am_nbead, &
soft_lambda,soft_alpha,soft_expo, &
AMOEBA_read_soft, soft_atom_range1,soft_atom_range2, soft_line, &
vdw_longrange_lambda
contains
soft_alpha,soft_lambda,soft_expo,vdw_longrange_lambda see descrition above
soft_atom_range1 and soft_atom_range2 are two arrays store the ranges of softcore atoms taken from soft_atom.txt
soft_line is the number of lines of soft_atom.txt
- AMOEBA_read_soft is a subroutine which read the soft_atom.txt file which will be called in amoeba_vdw.f.
- Include the new parameters in the amoeba_mdin namelist
subroutine AMOEBA_read_mdin(nf)
...skipping...
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,&
#ifndef DISABLE_AMOEBA_CG
dipole_scf_use_cg, dipole_scf_cg_niter, &
#endif /* DISABLE_AMOEBA_CG */
ee_dsum_cut,ee_damped_cut, &
thole_expon_coeff,vdw_taper,do_vdw_taper,do_vdw_longrange, &
compress,am_nbead,&
soft_lambda,soft_alpha,soft_expo,vdw_longrange_lambda
- Add subroutine AMOEBA_read_soft to read softcore atom list file
subroutine AMOEBA_read_soft()
integer pos_dash,length,i_range,i
character(len=20) atm_range
character(len=6) temp
!flag whether soft_atm.txt exists
logical alive
inquire(file='soft_atm.txt',exist=alive)
if(alive) then
do i=1,20
soft_atom_range1(i)=0;
soft_atom_range2(i)=0;
end do
soft_line=0
open (11,FILE='soft_atm.txt')
do while (.true.)
read(11,*,end=99) atm_range
!write(*,*) atm_range
pos_dash=scan(atm_range,'-')
!write(*,*) "position",pos_dash
length=len_trim(atm_range)
if (length > 0) then
soft_line=soft_line+1
endif
!write(*,*) "length",length
if(pos_dash.gt.0) then
temp=atm_range(1:pos_dash-1)
read(temp,*) soft_atom_range1(soft_line)
temp=atm_range(pos_dash+1:length)
read(temp,*) soft_atom_range2(soft_line)
else
read(atm_range,*) soft_atom_range1(soft_line)
read(atm_range,*) soft_atom_range2(soft_line)
endif
end do
else
write(6,*) 'soft_atm.txt not found'
call mexit(6,1)
endif
99 continue
close(11)
end subroutine AMOEBA_read_soft
amoeba_vdw.f
- softcore modification of vdw interaction
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, &
soft_lambda,soft_alpha,soft_expo, &
soft_atom_range1,soft_atom_range2, soft_line
integer,intent(in) :: i,ipairs(*),numtot
_REAL_,intent(in) :: xk,yk,zk,crd(3,*)
_REAL_,intent(inout) :: ene_vdw,frc(3,*),virial(3,3)
integer :: itran,it,jt,ih,jh,j,m,np,mask27,i_range
...skipping...
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
do i_range = 1,soft_line
if ((i.ge.soft_atom_range1(i_range).and. i.le.soft_atom_range2(i_range)) &
.xor. (j.ge.soft_atom_range1(i_range) .and.j.le.soft_atom_range2(i_range)) )then
eps = eps * soft_lambda ** soft_expo
t1 = (one + vdw_buf_delta)**7 / (soft_alpha * (1 - soft_lambda)**2 + (rho + vdw_buf_delta)**7)
t2 = (one + vdw_buf_gamma) / (soft_alpha * (1 - soft_lambda)**2 + rho7 + vdw_buf_gamma)
dt1drho = (-seven * (rho + vdw_buf_delta)**6 * t1) / (soft_alpha * (1 - soft_lambda)**2 &
+ (rho + vdw_buf_delta)**7)
dt2drho = (-seven * rho6 * t2) / (soft_alpha * (1 - soft_lambda)**2 + rho7 + vdw_buf_gamma)
endif
end do
f = eps*t1*(t2 - two)
dfdr = eps*(dt1drho*(t2 - two) + t1*dt2drho)*drhodr
- longrange correction softcore
subroutine AM_VDW_longrange_factor(num_atoms)
use amoeba_mdin, only : do_vdw_taper, &
soft_lambda,soft_atom_range1,soft_atom_range2,&
soft_line,vdw_longrange_lambda,AMOEBA_read_soft
use constants, only : zero,one,two,three,four,five,seven,pi
integer,intent(in) :: num_atoms
integer ier,n,nt,kdel,ndel,i,j,i_range,i_soft,i_soft_type
integer,save,allocatable :: lig_type_ct(:)
_REAL_ :: r,r1,r2,r3,r4,r5,req,eps,f,sume,sumv,t1,t2,rho,switch,delr
_REAL_ :: dt1drho,dt2drho,drhodr,dfdr,dswitch_dr,g1,g2
allocate(vdw_type_count(num_vdw_atom_types),stat=ier)
REQUIRE(ier==0)
allocate(lig_type_ct(num_vdw_atom_types),stat=ier)
REQUIRE(ier==0)
if (vdw_longrange_lambda .ne. 1.0 .or. soft_lambda .ne. 1.0) then
call AMOEBA_read_soft()
endif
do n = 1,num_vdw_atom_types
vdw_type_count(n) = 0
lig_type_ct(n) = 0
enddo
do n = 1,num_atoms
nt = vdw_atom_type(n)
vdw_type_count(nt) = vdw_type_count(nt) + 1
do i_range = 1, soft_line
if (n .ge. soft_atom_range1(i_range).and. n.le.soft_atom_range2(i_range)) then
lig_type_ct(nt) = lig_type_ct(nt) + 1
endif
enddo
enddo
...skipping...
! note the 2*pi below not 4*pi---since we do each i,j pair 2x
ene_vdw_longrange_factor = ene_vdw_longrange_factor + two*pi*sume* (vdw_type_count(i)*vdw_type_count(j) &
-(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j)))
vir_vdw_longrange_factor = vir_vdw_longrange_factor + two*pi*sumv*(vdw_type_count(i)*vdw_type_count(j) &
-(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j)))
enddo
enddo
end subroutine AM_VDW_longrange_factor
Modifying PMEMD 10
mdin_amoeba_dat.fpp
- Add the parameter(s) to the variable definitions in amoeba_mdin.f module mdin_amoeba_dat_mod and assign initial values to them:
module mdin_amoeba_dat_mod
#ifdef AMOEBA
implicit none
! Data that should be broadcast to slave processes from the master:
integer, parameter :: mdin_amoeba_int_cnt = 24
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
integer :: soft_atom_range1(20), soft_atom_range2(20), soft_line
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
common soft_atom_range1, soft_atom_range2, soft_line
save :: / mdin_amoeba_int /
integer, parameter :: mdin_amoeba_dbl_cnt = 10
double precision compress, dipole_scf_tol, ee_dsum_cut, &
ee_damped_cut, sor_coefficient, &
thole_expon_coeff, vdw_taper, &
soft_lambda, soft_alpha, vdw_longrange_lambda
common / mdin_amoeba_dbl / compress, dipole_scf_tol, ee_dsum_cut, &
ee_damped_cut, sor_coefficient, &
thole_expon_coeff, vdw_taper, &
soft_lambda, soft_alpha, vdw_longrange_lambda
save :: / mdin_amoeba_dbl /
contains
- assign the initial values for these softcore parameters in subroutine init_mdin_amoeba_dat
subroutine init_mdin_amoeba_dat
use file_io_dat_mod
use file_io_mod
implicit none
! Local variables:
integer :: ifind
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, &
soft_lambda, soft_alpha, soft_expo, vdw_longrange_lambda
...skipping...
soft_lambda=1.d0
soft_alpha=0.5d0
soft_expo=4
vdw_longrange_lambda=1.d0
- Add subroutine AMOEBA_read_soft
subroutine AMOEBA_read_soft
use parallel_dat_mod
integer pos_dash,length,i_range,i
character(len=20) atm_range
character(len=6) temp
logical alive
inquire(file='soft_atm.txt',exist=alive)
if (alive) then
do i=1,20
soft_atom_range1(i)=0;
soft_atom_range2(i)=0;
end do
soft_line=0
open (11,FILE='soft_atm.txt')
do while (.true.)
soft_line=soft_line+1
read(11,*,end=99) atm_range
pos_dash=scan(atm_range,'-')
length=len_trim(atm_range)
if(pos_dash.gt.0) then
temp=atm_range(1:pos_dash-1)
read(temp,*) soft_atom_range1(soft_line)
temp=atm_range(pos_dash+1:length)
read(temp,*) soft_atom_range2(soft_line)
else
read(atm_range,*) soft_atom_range1(soft_line)
read(atm_range,*) soft_atom_range2(soft_line)
endif
end do
else
write(6,*)'soft_atm.txt not found'
call mexit(6, 1)
endif
99 continue
close(11)
end subroutine AMOEBA_read_soft
amoeba_direct.fpp
softcore modification for vdw interactions in subroutine am_vdw_direct_ene_frc_i
subroutine am_vdw_direct_ene_frc_i(atm_i, img, ipairs_sublst, x_tran, &
pair_cnt, crd, ene_vdw, frc, img_frc, &
virial, img_atm_map)
use mdin_amoeba_dat_mod, only : do_vdw_taper, soft_lambda, soft_expo, soft_alpha, dipole_scf_iter_max, &
soft_atom_range1,soft_atom_range2, soft_line
use amoeba_flags_mod
use img_mod
implicit none
skipping
! Local variables:
integer :: itran, it, jt, ih, jh, idx
integer :: sublst_idx
integer :: atm_j, img_j
integer :: i_range
double precision :: wi, wj
...skipping...
if (delr2 .lt. vdw_switch_off_2) then
jt = vdw_atom_type(atm_j)
eps = vdw_epsilon(jt, it)
rad = vdw_radius(jt, it)
delr = sqrt(delr2)
rho = delr / rad
rho6 = rho**6
rho7 = rho6 * rho
t1 = ((1.d0 + vdw_buf_delta) / (rho + vdw_buf_delta))**7
t2 = (1.d0 + vdw_buf_gamma) / (rho7 + vdw_buf_gamma)
dt1drho = -7.d0 * t1 / (rho + vdw_buf_delta)
dt2drho = -7.d0 * t2 * (rho6 / (rho7 + vdw_buf_gamma))
drhodr = 1.d0 / rad
do i_range = 1,soft_line
if(((atm_i.ge.soft_atom_range1(i_range)).and.(atm_i.le.soft_atom_range2(i_range))) &
.or. ((atm_j.ge.soft_atom_range1(i_range)) .and.(atm_j.le.soft_atom_range2(i_range)))) then
eps = eps * soft_lambda ** soft_expo
t1 = (1.d0 + vdw_buf_delta)**7 / (soft_alpha * (1.d0 - soft_lambda)**2 + (rho + vdw_buf_delta)**7)
t2 = (1.d0 + vdw_buf_gamma) / (soft_alpha * (1.d0 - soft_lambda)**2 + rho**7 + vdw_buf_gamma)
dt1drho = (-7.d0 * (rho + vdw_buf_delta)**6 * t1) / (soft_alpha *(1.d0 - soft_lambda)**2 + (rho + vdw_buf_delta)**7)
dt2drho = (-7.d0 * rho**6 * t2) / (soft_alpha * (1.d0 - soft_lambda)**2 + rho**7 + vdw_buf_gamma)
endif
end do
f = eps * t1 * (t2 - 2.d0)
dfdr = eps * (dt1drho * (t2 - 2.d0) + t1 * dt2drho) * drhodr
The modification to the long-range correction is to separate out the ligand contribution:
total vdw type count =Other+Ligand (O+L)*(O+L)=O*O+O*L+L*O+L*L O=T-L, so O*L+L*O+L*L=(T-L)*L+L*(T*L)+L*L=T*L+L*T-L*L
amoeba_vdw.f
- modify subroutine am_vdw_longrange_factor
subroutine am_vdw_longrange_factor(atm_cnt)
use mdin_amoeba_dat_mod, only : do_vdw_taper, soft_lamda,AMOEBA_read_soft, &
soft_atom_range1,soft_atom_range2,soft_line,&
vdw_longrange_lambda
use gbl_constants_mod, only : PI
implicit none
! Formal arguments:
integer, intent(in) :: atm_cnt
! Local variables:
integer :: ier, n, nt, kdel, ndel, i, j
integer :: i_range,i_soft,i_soft_type
double precision :: r, r1, r2, r3, r4, r5 double precision :: req, eps, f, sume, sumv, t1, t2, rho, switch, delr double precision :: dt1drho, dt2drho, drhodr, dfdr, dswitch_dr, g1, g2
integer,save,allocatable :: lig_type_ct(:)
if (vdw_longrange_lambda .ne. 1.0 .or. soft_lamda .ne. 1.0) then
call AMOEBA_read_soft()
endif
allocate(lig_type_ct(vdw_param_cnt),stat=ier)
do n = 1, vdw_param_cnt
vdw_type_count(n) = 0
lig_type_ct(n) = 0
end do
do n = 1, atm_cnt
nt = vdw_atom_type(n)
vdw_type_count(nt) = vdw_type_count(nt) + 1
do i_range = 1, soft_line
if (n .ge. soft_atom_range1(i_range).and. n.le.soft_atom_range2(i_range)) then
lig_type_ct(nt) = lig_type_ct(nt) + 1
endif
enddo
end do
...skipping...
! Note the 2 * pi below not 4 * pi---since we do each i, j pair 2x.
ene_vdw_longrange_factor = ene_vdw_longrange_factor + 2.d0 * PI * sume * (vdw_type_count(i)*vdw_type_count(j) &
-(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j)))
vir_vdw_longrange_factor = vir_vdw_longrange_factor + 2.d0 * PI * sumv * (vdw_type_count(i)*vdw_type_count(j) &
-(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j)))
end do
end do
return
end subroutine am_vdw_longrange_factor
The modification to the long-range correction is to separate out the ligand contribution:
total vdw type count =Other+Ligand (O+L)*(O+L)=O*O+O*L+L*O+L*L O=T-L, so O*L+L*O+L*L=(T-L)*L+L*(T*L)+L*L=T*L+L*T-L*L
AMBER modification to use special damping factor of metal ions (Ca2+, Mg2+ and Zn2+)
SANDER 10
In Sander 10: There are three files modified: amoeba_interface.f, amoeba_induced.f and sander.f.
In tinker, we did: polarity(i)**sixth*pgmion(i) where pgmion is 1.6 for Mg++;
in amber, sq_polinv(i) is the sqrt(alpha_i); where alpha_i is the atomic polarizability
exp{-0.39*[r/(alphai^1/6*alphaj^1/6)]^3}
- 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,numpolar
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) - 24.305) < 0.0001d0 ) &
write (6,*) "Mg 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 ) &
sq_polinv(n) = 1 /1.60d0**3*sq_polinv(n)
if ( abs(mass(n) - 65.409) < 0.0001d0 ) &
write (6,*) "Zn damp",n,mass(n),polarizability(n)
if ( abs(mass(n) - 65.409) < 0.0001d0 ) &
sq_polinv(n) = 1 /1.23d0**3*sq_polinv(n)
endif
enddo
AM_INDUCED_readparm = 1
do_flag = ibset(do_flag,VALID_BIT)
call cg_init(num_atoms)
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
PMEMD 10
- In amoeba_induced.fpp, add
! The data previously in the sq_polinv array is now stored in atm_qterm()
! over in prmtop_dat_mod. This is the easiest way to switch between using
! a field in the img type for something useful in both amber pme and amoeba.
do n = 1, polar_atm_cnt
if (polarizability(n) .gt. polmin) then
is_polarizable(n) = .true.
atm_qterm(n) = 1.d0 / sqrt(polarizability(n))
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)
if ( abs(atm_mass(n) - 65.409) < 0.0001d0 ) &
write (6,*) "Zn damp",n,atm_mass(n),polarizability(n)
if ( abs(atm_mass(n) - 65.409) < 0.0001d0 ) &
atm_qterm(n) = 1 /1.23d0**3*atm_qterm(n)
else
is_polarizable(n) = .false.
atm_qterm(n) = 1.d0 / sqrt(polmin)
end if
ind_dip_d(:, n) = 0.d0
ind_dip_p(:, n) = 0.d0
end do
HFE Automation
How to use tinker and amber together (translated)
-added by xiaojia
most importantly, make sure you can get consistent energy from both software,
Correction "do_vdw_longrange=1" works if you add "VDW-CORRECTION" to the *.key for tinker
1) cut offs in Vdw and ewald for both should be the same.
2) "pme-grid" in Tinker.key (I don't have this in my mono.key...)= "nfft1=36,nfft2=36,nfft3=36," in amber.mdin
3) "polar-eps(default =10^-6)"in Tinker.key = "dipole_scf_tol" in amber.mdin
4) "polar-sor" in Tinker.key = "sor_coefficient" in amber.mdin
5) change "do_vdw_longrange=1" to "=0", because in thinker it is off
6) "ewald-alpha" in Tinker.key = "ew_coeff" in amber.mdin
7) "pme-order" in Tinker.key = the "order" in amber.mdin
Introduction
hfe_setup.pl
- This script is used to set up all free energy perturbation steps automatically to calculate hydration free energy.
- 1. Scale the electrostatic and VDW parameters first and then in TINKER4 and then convert all the related files in each step to amber format.
- 2. Setup all BAR analysis files for each step
- 3. Setup gas phase intramolecular recharging steps for the ligand
- 4. Equilibration of the lig-wat system for 100ps before starting FEP
Files needed:
a. xyz file of the small molecule (lig.xyz)
b. parameter file in TINKER4 format (amoeba.prm)
c. xyz file of the water box (watbox.xyz)
d. amber mdin template (tmp_amber.mdin) Note this format only works for sander 8 or 9. 10 uses iamoeba instead of do_amoeba. verbose is replaced by amoeba_verbose=0
short md, nve ensemble &cntrl irest=0, nstlim=500000, ntpr=500, ntwr=500,ntave=500, nscm=500,ntwx=500, dt=0.001, vlimit=10.0, cut=12.,maxcyc=50,ntmin=2,imin=0, ntt=1, temp0=298.0,tempi=298.0,tautp=0.5, / &ewald nfft1=36,nfft2=36,nfft3=36, 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.01,dipole_scf_iter_max=30, sor_coefficient=0.7,ee_damped_cut=4.5,ee_dsum_cut=6.8, /
e. optional: input file
lig.xyz #ligand xyz file amoeba.prm #tinker4 parameter file watbox.xyz #water box file 28.77 #size of the water box (a-axis) 12.0 #cutoff for simulation n #using default electrostatic scaling protocol? 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 # manuallt assign protocol n #using default vdw scaling protocol? 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 # manuallt assign protocol tmp_amber.mdin #amber mdin template
- Usage
$ hfe_setup.pl Enter Cartesian Coordinate File Name of the ligand: mes.xyz Input the parameter files(.prm): amoeba_2009_ver4.prm.isoprop Enter the file name of waterbox(.xyz): box30.xyz Input the size of the box: 28.77 Enter the cutoff you want: 12.0 Would you like to use the default electrostatic scaling scheme 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 [Y/N]? y Would you like to use the default Van der Waals scaling scheme 0.9 0.8 0.75 0.7 0.65 0.6 0.5 0.4 0.2 0 [Y/N]? n How do you want to scale the van der Waals(0.9 0.8 ... 0.0: 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Input the mdin file template: tmp_amber.mdin
For comvenience, you can make an input file (see above) and type
$hfe_setup.pl < input &
When you see
The equilibration takes approximately 2 days. Please start jobs until then. 27580
It means equilibration has already started. You can start simulation jobs for each step when equilibration is complete.
BAR-amber
- This program is for analyzing the free energy change of perturbation of a MD simulation using Bennet Acceptance Ratio. The program requires two energy output files for two frame sets from AMBER, totally four.
- Usage
BAR-amber ii.e ij.e ji.e jj.e start_frame end_frame guessed_init
ii.e: energy file with first trajectory and first prm file
ij.e: energy file with first trajectory and second prm file
ji.e: energy file with second trajectory and first prm file
jj.e: energy file with second trajectory and second prm file
BAR-tinker
- This program is for analyzing the free energy change of perturbation of gas phase intramolecular recharging process using Bennet Acceptance Ratio. The program requires two energy output files for two frame sets from TINKER, totally four.
- Usage
BAR-tinker ii.e ij.e ji.e jj.e start_frame end_frame guessed_init
ii.e: energy file with first trajectory and first prm file
ij.e: energy file with first trajectory and second prm file
ji.e: energy file with second trajectory and first prm file
jj.e: energy file with second trajectory and second prm file
amoeba_conv.pl
- This program is to covert parameter files between tinker4 and tinker5.
- Currently, AMBER only recognizes TINKER4 parameter format. So if you have a TINKER5 format parameter file, you need to convert it to TINKER4.
- Usage
amoeba_conv.pl input_param out_param 0/1(0 for tinker4-tinker5,1 for tinker5-tinker4)
Note
When we scale the electrostatic parameters and polarizability of solute to 0, we actually can not put exact 0 for all charge, dipoles and quadrupoles for that type in the prm file. Tinker would think sites of this type can be deleted (instead of having zero parameters) which cause problem when we do perturbation later. So we get around this by give a charge of 0.00001 (all other dipole and quadrupole are zeros). This way TINKER will keep it as a multipole site. But when you run "analyze * PC", it only prints out the first few zeros (this file is then converted into amber prmtop) so effectively this site has 0 electrostatics.
Examples and Test
Example files are in: http://biomol.bme.utexas.edu/~yues/download/example
Download
hfe_setup.pl:http://biomol.bme.utexas.edu/~yues/download/hfe_setup.pl
amoeba_conv.pl:http://biomol.bme.utexas.edu/~yues/download/amoeba_conv.pl
BAR-amber:http://biomol.bme.utexas.edu/~yues/download/BAR-amber
BAR-tinker:http://biomol.bme.utexas.edu/~yues/download/BAR-tinker
AMBER modification to use special damping factor of metal ions (Ca2+, Mg2+ and Zn2+)
SANDER
In Sander 9: There are three files modified: amoeba_interface.f, amoeba_induced.f and sander.f.
In tinker, we did: polarity(i)**sixth*pgmion(i) where pgmion is 1.6 for Mg++;
in amber, sq_polinv(i) is the sqrt(alpha_i); where alpha_i is the atomic polarizability
exp{-0.39*[r/(alphai^1/6*alphaj^1/6)]^3}
- 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,numpolar
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) - 24.305) < 0.0001d0 ) &
write (6,*) "Mg 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 ) &
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
iamoeba above is for sander 9; do_amoeba for sander 8
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
?
Softcore modifications to fixed-charge models in AMBER11
REMD modifications to AMOEBA in AMBER10
Tips running AMOBEA in AMBER
- When go from tinker and amber, tinker.xyz must match tinker.pdb (order etc.)
Use xyzpdb to generate a PDB that is consistent with xyz. But you have to do it twice, one with tinker.seq (which is produced early when pdbxyz); and the other is without the presence of tinker.seq. the former has correct residue id and name, but the atom order may be different than that xyz file. the later has correct order but incorrect residue name or ID.
Then combine the two using ~yues/bin/pdb_crct.pl
- Octahedron (when you use amoeba_parm generate inpcrd and prmtop from tinker, and you want to use octahedron:)
Change the last two lines in the generated inpcrd file.
Angle: from 90 to 109.4712; 0.109471219000E+03 0.109471219000E+03 0.109471219000E+03
In the gnerated prmtop file, two places need to be changed. In the last second row of FLAG POINTERS,
change the value 1 (the last third) to 2 meaning octahedron. %FLAG POINTERS %FORMAT(10I8)
1170 1 1 1 1 1 1 1 0 0 1170 390 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 2 0 0 0 0
%FLAG ATOM_NAME
ALso change the angle (search BOX) from 90 to 109 0.90E2 to 0.10947122E+03
- typical set up in tinker/pmemd/sander MD simulations
Example mdin files for ppmemd and sander: http://biomol.bme.utexas.edu/wiki/index.php/Research:Amber_tips_pr#direcories_in_.2Fhome.2Fpren.2FAMOEBA.2Fsolvation.2Fcompare.2F
vdw 12 (or 9) + LRC nfft = 1.2 * box size (kewald.f has the list, has ot be the product of powers of 2,3,5) polar-eps 0.01 polar-sor >= 0.6 (0.7 is good) do_longrange=1 tautp=0.5
PMEMD benchmark
Ranger and lonestar (infini band network) becnhamrk: 32-48 CPUs for large systems (50k atoms)
Gigabit cluster benchmark: Using 2x4-core compute nodes (used to be stampede on tacc)
pmemd v10 mpich2-fc7 23000 atoms, 1000 steps 1ps (Jenny)
corexnode cpu_time wall_time 8x1 0.35 hr 0.37 hr 8x2 0.21 hr 0.5 hr 8x4 0.17 hr 1.24 hr Lonestar16 0.15hr 0.15 hr Lonestar32 0.11 hr 0.14 hr Ranger16 0.20 hr 0.20 hr Ranger32 0.15 hr 0.15 hr
pmemd v10 mpich2-fc7 14000 atoms, 1000 steps 1ps (Jenny using yue's system)
oldnode4x1 0.205hr 0.207 hr newnode4x1 0.24 hr 0.24 hr newnode8x1 0.18 hr 0.18 hr
pmemd v10, mpich2-1.2.1p1-2.fc13, 10000 atoms, 1000 steps (1 ps)
cores CPU time wall time 4x1 0.17 hr 0.18 hr 8x1 0.14 hr 0.16 hr 4x2 0.12 hr 0.23 hr 8x2 0.11 hr 0.25 hr
Install pmemd.amoeba with Amber12
on bme-nova or water?
1. get the latest version of amber12 from /home/pren/amber12/yues/src/pmemd.amoeba/src or /home/yues/program/amber12/src/pmemd.amoeba/src
OLD: ~pren/amber-git/
2. set new AMBERHOME environment variable:
export AMBERHOME=/users/jzhen/amber12 # (bash) setenv AMBERHOME /users/jzhen/amber12 # (csh, tcsh) also add $AMBERHOME/bin to your PATH.
3. add Intel compiler and mpi path:
source /opt/intel/cce/10.0.023/bin/iccvars.csh source /opt/intel/fce/10.0.023/bin/ifortvars.csh setenv MPI_HOME /opt/mpich2-fc13-p26 setenv PATH $PATH':/opt/mpich2-fc13-p26/bin' setenv LD_LIBRARY_PATH /opt/intel/fce/10.0.023/lib:$LD_LIBRARY_PATH setenv LD_LIBRARY_PATH /opt/intel/cce/10.0.023/lib:$LD_LIBRARY_PATH source /opt/intel/fce/10.0.023/lib source /opt/intel/cce/10.0.023/lib
4. set C and Fortran compiler to intel:
setenv I_MPI_CC icc setenv I_MPI_CXX icpc setenv I_MPI_F77 ifort setenv I_MPI_F90 ifort
5. compile serial version first
cd $AMBERHOME ./configure intel cd src/pmemd.amoeba/src/ make -j 8 install
6. then compile mpi version
cd $AMBERHOME
./configure -mpi intel
or
./configure -mpi -nofftw3 intel (if you have trouble with fftw3 during configure). change config.h as below
important: the config.h use mpicc and mpif90, which is point to gcc and gfortran by default
please modify config.h in src/config.h (example /home/yues/program/amber12/src/config.h.mpi)
in #PMEMD Specific build flags,
add -f90=ifort to line PMEMD_F90=mpif90 -DMPI src/pmemd.amoeba/src/make install
change PMEMD_CC=mpicc to PMEMD_CC=icc
change CC=mpicc to CC=icc
Then
src/pmemd.amoeba/src/make install
7. add setenv LD_LIBRARY_PATH /opt/mpich2-fc13-p26/lib:$LD_LIBRARY_PATH
Install Amber10 pmemd.amba on Stampede
- The source code was taken from /home/yues/program/amber10/src/pmemd.amba/src-remd-respa
- Used the default module loaded, i.e. intel/13 and mvapich2/1.9
- The source on Stampede is located at
/work/01995/qw939/software/amber10/src/pmemd.amba/src-remd-respa
- The executable is located at
/work/01995/qw939/software/amber10/exe/pmemd.amba
1. Add following lines in configure file
32d31 < linux_stampede (Linux, Intel, x86_64) 394,397d392 < elif [ $platform = linux_stampede ]; then < cat << EOD >> config.h < MATH_LIBS = -L$MKL_HOME/lib/intel64 -lpthread < EOD
2. Then go to
/work/01995/qw939/software/amber10/src/pmemd.amba/config_data
and create a template called "linux_stampede.ifort.mvapich" by copy contents from "linux_em64t.ifort"
Then run
./configure linux_stampede ifort mvapich
This should give you a config.h file
3. Manually add the following lines to the config.h file if they are missing:
MATH_DEFINES = -DMKL MATH_LIBS = -L/opt/apps/intel/13/composer_xe_2013.2.146/mkl/lib/intel64 -mkl=parallel -lpthread MPI_HOME = /opt/apps/intel13/mvapich2/1.9 MPI_DEFINES = -DMPI -DSLOW_NONBLOCKING_MPI MPI_INCLUDE = -I$(MPI_HOME)/include MPI_LIBDIR2 = /opt/apps/intel13/mvapich2/1.9/lib MPI_LIBS = -L$(MPI_LIBDIR2) -lmpich -lrt -lpthread
4. Comment "integer :: iargc_wrap" in subroutine get_cmdline() in get_cmdline.fpp
5. This should be able to get you the pmemd.amba
6. Tested with examples in:
/home/yues/testamber/ele1.0/respa/pmemd12-verlet-respa/
and the REMD example:
/home/yues/program/amber12/test/rem_respa_amoeba/
AUTOMATION
Shi Yue wrote an AMOEBA automation script for the Comp Chem group in Genentech recently, which includes the automation of parameter generation for the ligand using Poltype, preparing the protein for simulation, and setting up the simulation.
She put a copy on our cluster in case someone is interested in using it. Please check the following directory for the script, readme file and the example.
/home/yues/automation/new