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

  1. The name of an input file, like the one above
  2. The name of the output file
  3. An Amber prmtop-file
  4. 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