[gmx-users] Implicit solvation: nonpolar term?

[rdwducl]

Hello,

I am using GBSA/OBC implicit solvation with an ACE-type approximation for
the nonpolar term. Having looked at the GROMACS code, it would appear to be
using the following form for the nonpolar free energy:

sum_i 4*pi*tension*(R_i + R_s)^2 * (R_i/b_i)^6

where R_i and R_s are constants (vdw and probe radii), and b_i is the
conformation-dependent Born radius. The idea is that as solutes come
together, b_i increases and so the free energy decreases.

To test this effect, I have written a script that iteratively creates
configs with the 7 atoms (of type "CH2"):

(0,0,0)
(r,0,0)
(-r,0,0)
(0,r,0)
(0,-r,0)
(0,0,r)
(0,0,-r)

and then computes the GB nonpolar energy as a function of r in [0.1,1.0] nm.
What I'm finding is a /constant/ nonpolar energy, i.e. it is not varying
with r!

My input files are below. Am I doing something wrong?

Kind regards.


pmf.top ----------------------------------

#include "ffG45a3.itp"

[ implicit_genborn_params ]
CH2 1.0 1 1.0 0.250 1.0

[ moleculetype ]
; molname       nrexcl
CH2             1

[ atoms ]
; id    at type         res nr  residu name     at name  cg nr  charge   
mass
1       CH2             1       CH2             CH2      1       0.00000 
14.000

[ system ]
implicit test case

[ molecules ]
CH2        7


pmf.mdp ----------------------------------

constraints         =  all-bonds
constraint_algorithm=  LINCS
integrator          =  md
dt                  =  0.00001
nsteps              =  1
nstxout             =  1
nstvout             =  1
nstfout             =  1
nstlog              =  1
nstenergy           =  1
nstxtcout           =  1
nstlist             =  1
ns_type             =  simple
coulombtype         =  cut-off
vdwtype             =  cut-off
rlist               =  0
rcoulomb            =  0
rvdw                =  0
pbc                 = no
optimize_fft        = yes
nstcomm             = 100
comm-grps           = system
; Berendsen temperature coupling is on 
tcoupl              =  no
; Energy monitoring
energygrps          =  CH2
; Semiisotropic pressure coupling is now on
Pcoupl              =  no
; Generate velocites is on at 300K.
gen_vel             =  no
; Implicit solvation setup
implicit_solvent    =  GBSA
gb_algorithm        =  OBC
nstgbradii          =  1
rgbradii            =  0
gb_epsilon_solvent  =  77.6
sa_algorithm        =  Ace-approximation
sa_surface_tension  = -1


run.sh --------------------------------------------------

#!/bin/bash

# handy function
function atom
{
    printf "%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" $1 CH2 CH2 $1 $2 $3 $4 >>
conf.gro
}

# loop from 0.1 nm to 1.0 nm
for((i=10; i<=100; i++))
do
    # r = i/100 nm
    r=$(echo "scale=10; 0.01*${i}"|bc)
    
    rm -f md.log

    # create the conf.gro file
    echo "test, t= 0.0" > conf.gro
    echo "7" >> conf.gro
    atom 1 $r 0.0 0.0
    atom 2 -$r 0.0 0.0
    atom 3 0.0 $r 0.0
    atom 4 0.0 -$r 0.0
    atom 5 0.0 0.0 $r
    atom 6 0.0 0.0 -$r
    atom 7 0.0 0.0 0.0
    printf "%10f%10f%10f\n" 10.0 10.0 10.0 >> conf.gro

    # run a quick single-point energy calculation
    grompp_mpi -f pmf.mdp -c conf.gro -p pmf.top > _E 2>&1
    mpirun -n 1 mdrun_mpi > _E 2>&1

    # compute and dump the nonpolar GB energy
    energy=$(grep -A 1 'GB Polarization' md.log|tail -n 1|awk '{print $2}')
    printf "%f\t%f\n" $r $energy
done

# cleanup
rm _E \#*



--
View this message in context: http://gromacs.5086.x6.nabble.com/Implicit-solvation-nonpolar-term-tp5013432.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.