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. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.