Continuum van der Waals

This document describes how to calculate continuum van der Waals energies with a modified version of the Amber pbsa program.

The program is callled cvdw and is located locally in /away/bio/AMBER/Amber11/AmberTools/src/cvdw/. The modification done are explained below.

Samuel Genheden, 2012.

Typical input-files

Here is a typical input file that will use 6/12 decomposition of Lennard-Jones potential, as well as SAV for the cavity term

6/12 decomposition, SAV for cavity term &cntrl imin=0,ntx=1,ipb=1,inp=2, / &pb cutres=12, fillratio=2, dprob=1.4, space=1.0, radiopt=1, decompopt=1, use_rmin=1, sprob=0.557, vprob=1.300, rhow_effect=1.129, use_sav=1, cavity_surften=0.0378, cavity_offset=-0.5692, /

It can be used for complex, receptor, and ligand.

If one would like to calculate the receptor energies but with a non-interacting ligand (P0), use the following input

6/12 decomposition, SAV for cavity term &cntrl imin=0,ntx=1,ipb=1,inp=2,zstart=7709,zstop=7720, / &pb cutres=12, fillratio=2, dprob=1.4, space=0.50, radiopt=1, decompopt=1, use_rmin=1, sprob=0.557, vprob=1.300, rhow_effect=1.129, use_sav=1, cavity_surften=0.0378, cavity_offset=-0.5692, /

and substitute the values of the zstart and zstop variables with the first and last atom number of the ligand.

Please, consult the manual of pbsa to find out the meaning of all the variables and what alternatives exists.

Execute cvdw

The program is run by typing

cvdw -O -i cvdw.in -o cvdw_com.out -p lig.prmtop -c lig.rst

where the arguments are

The name of an input file, like the one above
The name of the output file
An Amber prmtop-file
An Amber restart/coordinate-file

Output

The interesting output is found after the lines

-------------------------------------------------------------------------------- 4. RESULTS -------------------------------------------------------------------------------- Cavity solvation energy 26.0231 Dispersion solvation energy -48.3042 Repulsion solvation energy 20.9073

Note!: All energies are in kcal/mol.

Known limitations

A current limitation is that the repulsion energy will only be evaluated if one use decompopt=1, i.e. 6/12 decomposition. However, if one is only interested in the total van der Waals energy, this is no limitation.

I have not been able to successfully turn off Poisson-Boltzmann calculations completely. Therefore, errors might occur due to the instabilities of the PB-solver.

For the same reason, one might run into memory problems, and then one can try to run the calculations on a cluster.

Modifications to the pbsa-program

Modifications have been made to the np_force.f, pb_init.f, pb_read.f, md.h and Makefile files.

This is a diff between the modified and original np_force.f file.


33c33

< subroutine np_force( natom,nres,ntypes,ipres,iac,ico,cn1,cn2,x,f,enbrfcav,enbrfdis,enbrfrep )

---

> subroutine np_force(natom,nres,ntypes,ipres,iac,ico,cn1,cn2,x,f,enbrfcav,enbrfdis )

48d47

< #  integer zstart,zstop

56,61d54

< 

< #ifdef SANDER

< #  zstart = -1

< #  zstop = -1

< #endif

< 

66c59

<    _REAL_ enbrfcav, enbrfdis,enbrfrep

---

>    _REAL_ enbrfcav, enbrfdis

78c71

<    enbrfcav= ZERO; enbrfdis = ZERO ; enbrfrep = ZERO

---

>    enbrfcav= ZERO; enbrfdis = ZERO

86d78

<       write(6, '(1x,a,f12.4)') 'First atom radii=',mdsig(1)

123c115

<       do iatm = 1, natom      
     
---

>       do iatm = 1, natom

127,131c119

<             if ((zstart.gt.0).and.(iatm.ge.zstart).and.(iatm.le.zstop)) then

<               epsln_iatm = ZERO

<             else

<               epsln_iatm = cn2(ic)/(256.0d0*mdsig(iatm)**6)   

<             endif

---

>             epsln_iatm = cn2(ic)/(256.0d0*mdsig(iatm)**6)   

149d136

<          write(6, '(1x,a,f12.4)') 'Repulsion solvation energy', enbrfrep

200d186

<    ! _REAL_ enbrfrep          ! returned solvation repulsion energy

217c203

<    _REAL_ adis, arep

---

>    _REAL_ adis

222c208

<    _REAL_ tmpg, tmpf, gcorrec, tmpg2

---

>    _REAL_ tmpg, tmpf, gcorrec

263d248

<       if ((zstart.lt.0).and.(zstop.gt.0).and.((iatm.lt.abs(zstart)).or.(iatm.gt.zstop))) 
cycle

316d300

<          tmpg2 = ZERO

347d330

<                   arep =  a(iatm)/(THREE*((ris2*ris4)*(ris2*ris4)))

360d342

<                   arep = ZERO

373c355

<                   arep = ZERO

---

>          
       
377d358

<                tmpg2 = tmpg2 + arep*dotprot

415d395

<          enbrfrep = enbrfrep + tmpg2*delts(jatm)

This is a diff between the modified and the original pb_read.f file.


32c32

<                     idecomp,zstart,zstop

---

>                     idecomp

141,143d140

<    zstart = -1

<    zstop = -1

< 

295,298d291

< 

<    write(6,'(/a)') 'Dummy ligand definition:'

<    write(6,'(5x,2(a,i8))') 'zstart  =',zstart,&

<          ', zstop   =',zstop

832d824

<    ! SG 20/12 Make this possible

834,836c826,828

<       !write(6, *) 'PB Bomb in pb_read(): use of radi other than vdw sigma for'

<       !write(6, *) '                      np solvation dispersion/cavity is not supported!'

<       !call mexit(6, 1)

---

>       write(6, *) 'PB Bomb in pb_read(): use of radi other than vdw sigma for'

>       write(6, *) '                      np solvation dispersion/cavity is not supported!'

>       call mexit(6, 1)

839c831

<       !donpsa = .false.

---

>       donpsa = .false.

This is a diff between the modified and the original pb_init.f file.


330,343c330,337

<          ! SG 20/12 take radii from input instead

<          if (radiopt == 1) then

<            ic = ico(ntypes*(iac(iatm)-1) + iac(iatm))

<            if (cn2(ic) /= ZERO) then

<               mdsig(iatm) = (cn1(ic)/cn2(ic))**(SIXTH)/2 ! this is sigma

<               rmin(iatm) = mdsig(iatm)*(2.0d0**(SIXTH)) ! this is Rmin 

<            else

<               mdsig(iatm) = ZERO

<               rmin(iatm) = ZERO

<            endif

<         !else

<           !mdsig = rin(iatm)

<           !rmin = rin(iatm)

<         endif

---

>          ic = ico(ntypes*(iac(iatm)-1) + iac(iatm))

>          if (cn2(ic) /= ZERO) then

>             mdsig(iatm) = (cn1(ic)/cn2(ic))**(SIXTH)/2 ! this is sigma

>             rmin(iatm) = mdsig(iatm)*(2.0d0**(SIXTH)) ! this is Rmin 

>          else

>             mdsig(iatm) = ZERO

>             rmin(iatm) = ZERO

>          endif

This is a diff between the modified and the original md.h file.


21,22c21,23

<       ipb, inp, ibgion, ienion  !58

< parameter (BC_MDI=56)

---

>       ipb, inp, ibgion, ienion, & !58

>       zstart,zstop !60

> parameter (BC_MDI=58)

32c33

<       ipb,inp,ibgion,ienion

---

>       ipb,inp,ibgion,ienion,zstart,zstop

The Makefile was modified such that all strings of pbsa were replaced with strings of cvdw