>From ee048ea9abbb76d3b6e52fd372aae57ca7110557 Mon Sep 17 00:00:00 2001 From: wmsmith Date: Mon, 30 Sep 2013 13:48:50 -0700 Subject: [PATCH] Minor changes to pbsa and mmpbsa to fix minor bugs and make backward compatible with amber13 defaults. --- .../src/mmpbsa_py/MMPBSA_mods/createinput.py | 4 ++-- AmberTools/src/mmpbsa_py/mmpbsa_energy.nab | 2 +- AmberTools/src/pbsa/force.F90 | 2 +- AmberTools/src/pbsa/pb_init.F90 | 24 +++++++++++++++++--- AmberTools/src/pbsa/pb_read.F90 | 10 ++++---- .../test/mmpbsa_py/02_Poisson_Boltzmann/Run.PB | 2 +- 6 files changed, 31 insertions(+), 13 deletions(-) diff --git a/AmberTools/src/mmpbsa_py/MMPBSA_mods/createinput.py b/AmberTools/src/mmpbsa_py/MMPBSA_mods/createinput.py index 1ccae9b..45cfd99 100644 --- a/AmberTools/src/mmpbsa_py/MMPBSA_mods/createinput.py +++ b/AmberTools/src/mmpbsa_py/MMPBSA_mods/createinput.py @@ -240,7 +240,7 @@ class GBNabInput(NabInput): class PBNabInput(NabInput): """ The PB input file for the mmpbsa_py_energy NAB program """ - input_items = {'inp' : 2, 'smoothopt' : 1, 'radiopt' : 1, + input_items = {'inp' : 2, 'smoothopt' : 1, 'radiopt' : 0, 'npbopt' : 0, 'solvopt' : 1, 'maxitn' : 1000, 'nfocus' : 2, 'bcopt' : 5, 'eneopt' : 2, 'fscale' : 8, 'epsin' : 1.0, 'epsout' : 80.0, @@ -317,7 +317,7 @@ class SanderPBSAInput(SanderInput): 'imin':5, 'maxcyc':1, 'igb':10, 'ipb':2, 'inp':2, 'idecomp':0, 'dbfopt':1, 'epsin':1.0, 'epsout':80.0, - 'istrng':0.0, 'radiopt':1, 'sprob':0.557, 'dprob':1.4, + 'istrng':0.0, 'radiopt':0, 'sprob':0.557, 'dprob':1.4, 'space':0.5, 'maxitn':1000, 'cavity_surften':0.0378, 'cavity_offset':-0.5692, 'fillratio':4.0, 'ioutfm':0, 'dec_verbose':0} diff --git a/AmberTools/src/mmpbsa_py/mmpbsa_energy.nab b/AmberTools/src/mmpbsa_py/mmpbsa_energy.nab index 91cdf33..67ca4bc 100644 --- a/AmberTools/src/mmpbsa_py/mmpbsa_energy.nab +++ b/AmberTools/src/mmpbsa_py/mmpbsa_energy.nab @@ -156,7 +156,7 @@ int SetDefaults() { pbopt.runpb = 0; pbopt.inp = 2; pbopt.smoothopt = 1; - pbopt.radiopt = 1; + pbopt.radiopt = 0; pbopt.npbopt = 0; pbopt.solvopt = 1; pbopt.maxitn = 100; diff --git a/AmberTools/src/pbsa/force.F90 b/AmberTools/src/pbsa/force.F90 index 0ca0cde..a180c37 100644 --- a/AmberTools/src/pbsa/force.F90 +++ b/AmberTools/src/pbsa/force.F90 @@ -174,7 +174,7 @@ subroutine force(xx,ix,ih,ipairs,x,f,ener,vir,r_stack,i_stack, & call pbtimer_stop(PBTIME_PBFORCE) end if ! ( ipb >= 1 ) - if ( inp /= 0 ) then + if ( inp /= 0 .and. ipb >= 1 ) then call pbtimer_start(PBTIME_NPFORCE) esurf = 0.0d0; edisp = 0.0d0 call np_force(natom,nres,ntypes,ix(i02),ix(i04),ix(i06),& diff --git a/AmberTools/src/pbsa/pb_init.F90 b/AmberTools/src/pbsa/pb_init.F90 index d0b9ee6..beb2140 100644 --- a/AmberTools/src/pbsa/pb_init.F90 +++ b/AmberTools/src/pbsa/pb_init.F90 @@ -161,7 +161,7 @@ subroutine pb_init(ifcap,natom,nres,ntypes,nbonh,nbona,ipres,iac,ico,numex,natex integer ires, iatm, jatm, maxmax, ic, i, j, jp, idum integer alloc_err(64) _REAL_ maxnba_l - _REAL_ rinchk + _REAL_ rinchk, crgchk _REAL_ ucrgh(natom), ucrga(natom) character (len=4) :: residue, resid(natom) integer focusresidue(nres) @@ -347,7 +347,9 @@ subroutine pb_init(ifcap,natom,nres,ntypes,nbonh,nbona,ipres,iac,ico,numex,natex end do if (rinchk == ZERO) then write(6,*) 'PB Bomb in pb_init(): Requested radi to be read in, but found none' - call mexit(6,1) + write(6,*) 'PB Bomb in pb_init(): No PB calculation done' + !call mexit(6,1) + ipb = 0 end if radi = rin ! for pb mdsig = rin ! for np @@ -358,8 +360,24 @@ subroutine pb_init(ifcap,natom,nres,ntypes,nbonh,nbona,ipres,iac,ico,numex,natex call phi_aaradi( natom, isymbl, radi ) else write(6,*) 'PB Bomb in pb_init(): Unknown radi assigment option', radiopt - call mexit(6,1) + write(6,*) 'PB Bomb in pb_init(): No PB calculation done' + !call mexit(6,1) + ipb = 0 end if + + ! check if atomic charges are correct input from prmtop file + + crgchk = ZERO + do iatm = 1, natom + crgchk = crgchk+acrg(iatm) + end do + if (crgchk == ZERO) then + write(6,*) 'PB Bomb in pb_init(): Requested charges to be read in, but found none' + write(6,*) 'PB Bomb in pb_init(): No PB calculation done' + !call mexit(6,1) + ipb = 0 + end if + if ( pbverbose ) then write(6,*) ' no. of atoms processed in PB initialization:', natom write(6, '(a5,3a6,5a10)') 'NUM','RESI','NAME','TYPE',& diff --git a/AmberTools/src/pbsa/pb_read.F90 b/AmberTools/src/pbsa/pb_read.F90 index 139a0f5..3fb10cd 100644 --- a/AmberTools/src/pbsa/pb_read.F90 +++ b/AmberTools/src/pbsa/pb_read.F90 @@ -673,7 +673,7 @@ subroutine pb_read srsas = .true. smoothopt = 1 - radiopt = 1 + radiopt = 0 radiscale = 1.0d0 npopt = 2 decompopt = 2 @@ -1005,10 +1005,10 @@ subroutine pb_read ! write(6, *) 'PB Info in pb_read(): npopt has been overwritten with inp.' npopt = inp endif - if ( inp == 2 .and. radiopt /= 1 ) then - write(6, *) 'PB Warning in pb_read(): radiopt has been reset to 1 because inp=2.' - radiopt = 1 - end if +! if ( inp == 2 .and. radiopt /= 1 ) then +! write(6, *) 'PB Warning in pb_read(): radiopt has been reset to 1 because inp=2.' +! radiopt = 1 +! end if if ( inp == 1 .and. use_sav == 1 ) then write(6, *) 'PB Warning in pb_read(): the cavity energy should be correlated to & the SASA when inp=1, so use_sav is reset to 0.' diff --git a/AmberTools/test/mmpbsa_py/02_Poisson_Boltzmann/Run.PB b/AmberTools/test/mmpbsa_py/02_Poisson_Boltzmann/Run.PB index 90c0580..a421917 100755 --- a/AmberTools/test/mmpbsa_py/02_Poisson_Boltzmann/Run.PB +++ b/AmberTools/test/mmpbsa_py/02_Poisson_Boltzmann/Run.PB @@ -28,7 +28,7 @@ MMPBSA input file for running PB keep_files=0, debug_printlevel=2 / &pb - istrng=0.1, + istrng=0.1,radipt=1 / EOF -- 1.7.9.5