!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! PBSA ! ! An analysis program for solvent-mediated electrostatic and non-electrostatic ! interactions in biomolecules. ! ! Please acknowledge your use of PBSA by citing: ! ! Luo, David, and Gilson, J. Comp. Chem. 23:1244-1253, 2002. ! ! Major Developers: ! ! Jun Wang, linear Poisson-Boltzmann numerical solvers ! Qin Cai, nonlinear Poisson-Boltzmann numerical solvers ! Xiang Ye, electrostatic energy and force numerical algorithms ! Meng-Juei Hsieh, program interface and parallel implementations ! Chuck Tan, non-electrostatic energy and force numerical algorithms ! Ray Luo, coordinator of overall development ! ! Additional contributing authors are listed in the code documentation. ! ! Departments of Molecular Biology and Biochemistry and Biomedical Engineering, ! University of California, Irvine, California ! ! Copyright (c) 2004-2009. The Regents of the University of California. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Portion of the routine "pb_exmol" is modified from the UHBD program ! (Comp. Phys. Comm. 91:57-95, 1995), copyrighted by University of Houston, ! 1989-2009. See program documentation for more information. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Overview: ! ! PBSA models the electrostatic solvation interaction by the Poisson-Boltzmann ! equation. The implementation uses the finite-difference method to solve the ! partial differential equation. Both linear and full nonlinear numerical ! solvers are implemented. Please refer to the following publications: ! ! Luo, David, and Gilson, J. Comp. Chem. 23:1244-1253, 2002. ! Cai, Wang, Zhao, and Luo, J. Chem. Phys. 130:14501, 2009. ! Wang and Luo, J. Comp. Chem. In Press, 2010. ! ! for implementation details of the linear solvers. The full nonlinear solvers ! are documented in: ! ! Cai, Wang, and Luo, J. Chem. Comp. Theo. In Press, 2010. ! ! The electrostatic energy and forces are computed based on the finite- ! difference grid potentials as discussed in the following publications: ! ! Lu and Luo, J. Chem. Phys. 119:11035-11047, 2003. ! Xiang, Wang, and Luo, J. Phys. Chem. Submited. 2010. ! ! The dielectric models and molecular surfaces used in the electrostatic ! solvation model are documented in: ! ! Ye, Wang, and Luo, J. Chem. Theo. Comp. In Press, 2010. ! Wang and Luo, J. Chem. Theo. Comp. Submitted, 2009. ! ! The parameters used in the electrostatic solvation model are documented in ! the following publication: ! ! Tan, Yang and Luo, J. Phys. Chem. 110:18680-18687, 2006. ! ! PBSA models the non-electrostatic solvation interaction by two separate ! terms in this release: dispersion (or van der Waals or attractive) and cavity ! (or hydrophobic or repulsive). The dispersion term is computed by a numerical ! integration over the solvent accessible surface area. The cavity term is ! modeled by a term proportional to the molecular surface or volume with a ! single proportional constant. See the following publication for details: ! ! Tan, Tan and Luo, J. Phys. Chem. 111:12263-12274, 2007. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! This file is part of PBSA. ! ! PBSA is free software; you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation; either version 2 of the License, or (at your option) any later ! version. ! ! PBSA is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS ! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PBSA; if not, write to the Free Software Foundation, Inc., 59 Temple Place, ! Suite 330, Boston, MA 02111-1307, USA. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ PBMD module parsing and initialization subroutine pb_read(smoothopt_, radiopt_, npbopt_, solvopt_, maxitn_, & nbuffer_, nfocus_, fscale_, npbgrid_, dbfopt_, bcopt_, scalec_, & eneopt_, frcopt_, nsnbr_, phiout_, phiform_, npbverb_, npopt_, & decompopt_, use_rmin_, use_sav_, maxsph_, maxarc_, ndofd_, ndosas_, & mpopt_, lmax_, pbprint_, & epsin_, epsout_, istrng_, pbtemp_, dprob_, iprob_, & accept_, fillratio_, space_, arcres_, cutres_, cutfd_, cutnb_, & sprob_, vprob_, rhow_effect_, cavity_surften_, cavity_offset_, & cutsa_, fmiccg_, ivalence_, laccept_, wsor_, lwsor_, radinc_, & expthresh_, offx_, offy_, offz_, sepbuf_) ! Module variables use poisson_boltzmann use dispersion_cavity use solvent_accessibility implicit none ! Common variables !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+pb_def.h PBSA cpp macro-definitions ! maximum number of neighboring atoms per atom ! maximum levels of focusing in fdpb ! maximum number of grid charges for boundary conditions !+End of pb_def.h !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Specification and control of Amber's Input/Output ! File names character(len=512) groupbuffer character(len=80) mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe,redir(8), & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version character owrite common /files/ groupbuffer, mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe, owrite, & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version ! put this in a seperate common block to stop the compiler from ! complaining about misalignment integer numgroup common/nmgrp/ numgroup ! File units ! An I/O Unit resource manager does not exist. integer MDCRD_UNIT integer MDEN_UNIT integer MDINFO_UNIT integer MDVEL_UNIT parameter ( MDINFO_UNIT = 7 ) parameter ( MDCRD_UNIT = 12 ) parameter ( MDEN_UNIT = 15 ) parameter ( MDVEL_UNIT = 13 ) integer, parameter :: CNSTPH_UNIT = 18, CPOUT_UNIT = 19 ! 18 was picked because CNSTPH uses it; conflicts are not expected. !integer MMTSB_UNIT !parameter ( MMTSB_UNIT = 18 ) ! File related controls and options character(len=80) title,title1 common/runhed/ title, title1 logical mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version common/mdin_flags/mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version integer BC_HULP ! size in integers of common HULP parameter ( BC_HULP = 9 ) integer ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ,mdm_nstep & !#endif ; ! line terminator for free-form version common/hulp/ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ! this variable is last in the common and does not increment BC_HULP ! ,mdm_nstep & !#endif ; ! line terminator for free-form version ! NMRRDR : Contains information about input/output file redirection ! REDIR and IREDIR contain information regarding ! LISTIN, LISTOUT, READNMR, NOESY, SHIFTS, DUMPAVE, ! PCSHIFT and DIPOLE respectively. If IREDIR(I) > 0, ! then that input/output has been redirected. integer iredir(8) common/nmrrdr/redir,iredir !-------------BEGIN md.h ------------------------------------------------ integer BC_MDI ! size in integers of common block mdi integer BC_MDR ! size in _REAL_s of common block mdr ! ... integer variables: integer nrp,nspm,ig,ntx,ntcx, &!5 ntxo,ntt,ntp,ntr,init, &!10 ntcm,nscm,isolvp,nsolut,klambda, &!15 ntc,ntcc,ntf,ntid,ntn, &!20 ntnb,nsnb,ndfmin,nstlim,nrc, &!25 ntrx,npscal,imin,maxcyc,ncyc, &!30 ntmin,irest,jfastw, &!33 ibgwat,ienwat,iorwat, &!36 iwatpr,nsolw,igb,iyammp, &!40 gbsa,vrand,iwrap,nrespa,irespa,nrespai,icfe, &!47 rbornstat,ivcap,iconstreff, &!50 idecomp,icnstph,ntcnstph,maxdup, &!54 ipb, inp, ibgion, ienion !58 parameter (BC_MDI=56) common/mdi/nrp,nspm,ig, & ntx,ntcx,ntxo,ntt,ntp,ntr,init,ntcm,nscm, & isolvp,nsolut,ntc,ntcc,ntf,ntid,ntn,ntnb,nsnb,ndfmin, & nstlim,nrc,ntrx,npscal,imin,maxcyc,ncyc,ntmin, & irest,jfastw,ibgwat,ienwat,iorwat, & iwatpr,nsolw,igb,iyammp,gbsa,vrand, & iwrap,nrespa,irespa,nrespai,icfe,rbornstat, & ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup, & ipb,inp,ibgion,ienion ! ... floats: double precision t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & !9 tol,taur,dx0,drms,vlimit,rbtarg(8),tmass,tmassinv, & !24 kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & !31 cut_inner,clambda,saltcon, & !34 solvph,rgbmax,fsmax,restraint_wt !38 parameter (BC_MDR=38) common/mdr/t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & tol,taur,dx0,drms,vlimit,rbtarg,tmass,tmassinv, & kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & cut_inner,clambda,saltcon, & solvph,rgbmax,fsmax,restraint_wt ! ... strings: character(len=4) iwtnm,iowtnm,ihwtnm character(len=256) restraintmask,bellymask common/mds/ restraintmask,bellymask,iwtnm,iowtnm,ihwtnm(2) !-------------END md.h ------------------------------------------------ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ pb_md.h PB variables shared with calling routines logical pbverbose logical pbprint logical pbgrid logical pbinit common /pb_mdl/ pbverbose common /pb_mdl/ pbprint common /pb_mdl/ pbgrid common /pb_mdl/ pbinit integer npbstep integer nsaslag integer npbgrid integer nsnbr integer nsnba integer ntnbr integer ntnba integer dofd integer ndofd integer dosas integer ndosas common /pb_mdi/ npbstep common /pb_mdi/ nsaslag common /pb_mdi/ npbgrid common /pb_mdi/ nsnbr common /pb_mdi/ nsnba common /pb_mdi/ ntnbr common /pb_mdi/ ntnba common /pb_mdi/ dofd common /pb_mdi/ ndofd common /pb_mdi/ dosas common /pb_mdi/ ndosas !+ End of pb_md.h !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Local variables integer npbverb, l, phiout, scalec double precision space ! Passed variables integer smoothopt_, radiopt_, npbopt_, solvopt_, maxitn_, & nbuffer_, nfocus_, fscale_, npbgrid_, dbfopt_, bcopt_,& scalec_, eneopt_, frcopt_, nsnbr_, phiout_, phiform_, & npbverb_, npopt_, decompopt_, use_rmin_, use_sav_, & maxsph_, maxarc_, ndofd_, ndosas_, mpopt_, lmax_ logical pbprint_ double precision epsin_, epsout_, istrng_, pbtemp_, dprob_, iprob_, & accept_, fillratio_, space_, arcres_, cutres_, cutfd_,& cutnb_, sprob_, vprob_, rhow_effect_, cavity_surften_,& cavity_offset_, cutsa_, fmiccg_, ivalence_, laccept_, & wsor_, lwsor_, radinc_, expthresh_, offx_, offy_, & offz_, sepbuf_ ! initialize with default values outphi = .false. phiout = 0 phiform = 0 epsin = 1.0d0 epsout = 80.0d0 istrng = 0.0d0 ivalence = 1.0d0 pbtemp = 300.0d0 scalec = 0 scalerf = .false. srsas = .true. smoothopt = 0 radiopt = 1 npopt = 2 decompopt = 1 use_rmin = 0 use_sav = 0 rhow_effect = 1.0d0 maxsph = 400 maxarcdot= 0 maxarc = 256 nbuffer = 0 dprob = 0.00d0 sprob = 1.60d0 vprob = 1.28d0 iprob = 2.00d0 radinc = sprob*0.5d0 expthresh = 0.2d0 arcres = 0.0625d0 cavity_surften = 0.04356d0 cavity_offset = -1.008d0 nfocus = 2 fscale = 8 level = 1 bcopt = 5 space = 0.5d0 savxm(1) = 0 savym(1) = 0 savzm(1) = 0 savxmym(1) = 0 savxmymzm(1) = 0 savh(1) = 0 savbcopt(1) = 0 savxm(nfocus) = 0 savym(nfocus) = 0 savzm(nfocus) = 0 savxmym(nfocus) = 0 savxmymzm(nfocus) = 0 savh(nfocus) = 0 savbcopt(nfocus) = 0 offx = 0.0d0 offy = 0.0d0 offz = 0.0d0 fillratio= 2.0d0 npbopt = 0 solvopt = 1 fmiccg = -0.30d0 accept = 1.0d-3 laccept = 0.1d0 wsor = 1.9d0 lwsor = 1.95d0 maxitn = 100 pbverbose = .false. pbprint = .true. pbgrid = .true. pbinit = .true. npbverb = 0 npbgrid = 1 ndofd = 1 dofd = 1 donpsa = .true. ndosas = 1 nsaslag = 100 nsnbr = 1 nsnba = 1 ntnbr = 1 ntnba = 1 cutres = 99.0d0 cutfd = 5.0d0 cutnb = 0.0d0 cutsa = 9.0d0 lastp = 0 pbgamma_int = 1.0 pbgamma_ext = 65.0 sepbuf = 4.0d0 mpopt = 0 lmax = 80 dbfopt = -1 eneopt = -1 frcopt = 0 ! reading parameters smoothopt=smoothopt_ radiopt=radiopt_ npbopt=npbopt_ solvopt=solvopt_ maxitn=maxitn_ nbuffer=nbuffer_ nfocus=nfocus_ fscale=fscale_ npbgrid=npbgrid_ dbfopt=dbfopt_ bcopt=bcopt_ scalec=scalec_ eneopt=eneopt_ frcopt=frcopt_ nsnbr=nsnbr_ phiout=phiout_ phiform=phiform_ npbverb=npbverb_ npopt=npopt_ decompopt=decompopt_ use_rmin=use_rmin_ use_sav=use_sav_ maxsph=maxsph_ maxarc=maxarc_ ndofd=ndofd_ ndosas=ndosas_ mpopt=mpopt_ lmax=lmax_ pbprint=pbprint_ epsin=epsin_ epsout=epsout_ istrng=istrng_ pbtemp=pbtemp_ dprob=dprob_ iprob=iprob_ accept=accept_ fillratio=fillratio_ space=space_ arcres=arcres_ cutres=cutres_ cutfd=cutfd_ cutnb=cutnb_ sprob=sprob_ vprob=vprob_ rhow_effect=rhow_effect_ cavity_surften=cavity_surften_ cavity_offset=cavity_offset_ cutsa=cutsa_ fmiccg=fmiccg_ ivalence=ivalence_ laccept=laccept_ wsor=wsor_ lwsor=lwsor_ radinc=radinc_ expthresh=expthresh_ offx=offx_ offy=offy_ offz=offz_ sepbuf=sepbuf_ ! check and post processing of options if ( npbverb > 0 ) pbverbose = .true. dofd = ndofd dosas = ndosas epsin = epsin*eps0 epsout = epsout*eps0 istrng = fioni * istrng pbkappa = SQRT( 2.0d0 * istrng / (epsout * pbkb * pbtemp) ) ! set up solvent probes for different solvation components ! if end user requests different probes for different solvation ! components, no need to set up different probes here. if ( dprob == 0.0d0 ) dprob = sprob ! set up solvent accessible arc resolutions and limits if ( maxarcdot == 0) maxarcdot = 1200*nint(0.5d0/arcres) ! check dprob and iprob if ( dprob >= iprob ) then write(6, *) 'PB Bomb in pb_read(): solvent probe cannot be smaller than ion prob' call mexit(6, 1) end if if ( dprob <= space .and. ipb == 1 ) then write(6, *) 'PB Warning in pb_read(): setting grid spacing larger than solvent probe' write(6, *) 'may cause numerical instability if ipb=1' call mexit(6, 1) end if if ( dprob <= space .and. ipb == 2 ) then write(6, *) 'PB Bomb in pb_read(): solvent probe cannot be smaller than grid spacing if ipb=2' call mexit(6, 1) end if ! check pb options if ( igb == 10 .and. ipb /= 1 ) then write(6, *) 'PB Info in pb_read(): igb has been overwritten by ipb' endif ! check nonpolar options if ( inp /= npopt ) then write(6, *) 'PB Info in pb_read(): npopt has been overwritten with inp' npopt = inp endif if ( inp == 2 .and. radiopt == 0 ) then 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) end if if ( radiopt == 0 ) then donpsa = .false. end if ! check force options if ( eneopt == -1 ) then if ( dbfopt == 0 ) then eneopt = 1 else if ( dbfopt == 1 .or. dbfopt == -1 ) then eneopt = 2 else write(6, *) 'PB Bomb in pb_read(): only dbfopt= 0 or 1 are supported. dbfopt is replaced by eneopt' call mexit(6, 1) end if else if ( dbfopt == 0 .and. eneopt == 1 ) then write(6, *) 'PB Info in pb_read(): dbfopt has been overwritten by eneopt' else if ( dbfopt == 1 .and. eneopt == 2 ) then write(6, *) 'PB Info in pb_read(): dbfopt has been overwritten by eneopt' else if ( dbfopt /= -1 ) then write(6, *) 'PB Bomb in pb_read(): dbfopt is ignored' end if end if if ( eneopt < 1 .or. eneopt > 2 ) then write(6, *) 'PB Bomb in pb_read(): only eneopt= 1 or 2 are supported' call mexit(6, 1) end if if ( frcopt < 0 .or. eneopt > 3 ) then write(6, *) 'PB Bomb in pb_read(): only frcopt= 0-3 are supported' call mexit(6, 1) end if if ( eneopt == 1 .and. ( frcopt == 2 .or. frcopt == 3 ) ) then write(6, *) 'PB Bomb in pb_read(): this combination of eneopt and frcopt is not supported' call mexit(6, 1) end if if ( eneopt == 2 .and. frcopt == 1 ) then write(6, *) 'PB Bomb in pb_read(): this combination of eneopt and frcopt is not supported' call mexit(6, 1) end if if ( ( frcopt == 2 .or. frcopt == 3 ) .and. bcopt == 5 ) then bcopt = 6 write(6, *) 'PB Info in pb_read(): bcopt has been reset from 5 to 6 for frcopt= 2 or 3' end if if ( frcopt == 1 .and. bcopt == 6 ) then bcopt = 5 write(6, *) 'PB Info in pb_read(): bcopt has been reset from 6 to 5 for frcopt=1' end if if ( frcopt > 0 .and. smoothopt == 0 ) then smoothopt = 1 write(6, *) 'PB Info in pb_read(): smoothopt has been reset to 1 for force compuation' end if if ( ivcap /= 0 .and. (eneopt /= 1 .or. frcopt /= 0) ) then write(6,*) "cap water should only be run with eneopt=1 and frcopt=0" call mexit(6,1) endif ! check numerical solution options if ( npbopt == 0 .and. solvopt > 4 ) then write(6, *) 'PB Bomb in pb_read(): solvopt>4 cannot be used to solve linear PB equation' call mexit(6, 1) endif if ( solvopt == 2 .and. bcopt == 6 ) then write(6, *) 'PB Bomb in pb_read(): bcopt=6 cannot be used with solvopt=2' call mexit(6, 1) end if if ( npbopt == 1 .and. solvopt > 6 ) then write(6, *) 'PB Bomb in pb_read(): nonsupported solvopt (>6) for e nonlinear PB equation' call mexit(6, 1) endif if ( npbopt == 1 .and. eneopt == 2 ) then eneopt = 1 write(6, *) 'PB Info in pb_read(): eneopt has been reset to be 1 for nonlinear PB equation' end if if ( npbopt == 1 .and. frcopt > 0 ) then write(6, *) 'PB Bomb in pb_read(): force computatoin is not supported for nonlinear PB equation' call mexit(6, 1) end if ! check cutoff options if ( cutfd > cutsa ) then cutsa = cutfd end if if ( cutnb /= 0 .and. cutfd > cutnb ) then cutnb = cutfd write(6, *) 'PB Info in pb_read(): cutnb has been reset to be equal to cutfd' end if if ( cutnb /= 0 .and. bcopt == 6 ) then cutnb = 0 write(6, *) 'PB Info in pb_read(): cutnb has been reset to 0 with bcopt=6', cutnb, bcopt end if if ( max(cutnb,cutsa,cutfd) > cutres ) then write(6, *) 'PB Bomb in pb_read(): cutnb/cutfd must be <= cutres', cutnb, cutfd, cutres call mexit(6, 1) end if if ( cutnb == 0 .and. eneopt == 1 .and. bcopt /= 6 ) then write(6, *) 'PB Bomb in pb_read(): cutnb=0 cannot be used with eneopt=1', cutnb, eneopt call mexit(6, 1) end if cutfd = cutfd**2 cutsa = cutsa**2 cutnb = cutnb**2 cutres = cutres**2 ! set buffer zone between the fine FD grid boundary and the solute surface: if ( nbuffer == 0 ) then if ( istrng == 0.0d0 ) then nbuffer = int(2.0d0*sprob/space)+1 else if ( sprob >= iprob ) then nbuffer = int(2.0d0*sprob/space)+1 else nbuffer = int(2.0d0*iprob/space)+1 end if end if if ( nbuffer >= fscale ) then nbuffer = 2*nbuffer+1 else nbuffer = 2*fscale+1 end if end if ! set flag to scale induced surface charges: if ( scalec == 1) scalerf = .true. ! set saved grid options savbcopt(1) = bcopt do l = 2, nfocus savbcopt(l) = 0 end do savh(nfocus) = space do l = nfocus - 1, 1, -1 savh(l) = savh(l+1)*fscale end do do l = 1, nfocus savxm(l) = 1 savym(l) = 1 savzm(l) = 1 savxmym(l) = 1 savxmymzm(l) = 1 end do ! set phimap output options when requested if ( phiout == 1 ) then outphi = .true. radiopt = 2 donpsa = .false. npopt = 0 end if end subroutine pb_read subroutine myechoin (ilun,iout) ! Rewritten by: Meng-Juei Hsieh ! A rewrite of echoin implicit none integer ilun, iout logical boolinside, boolhastag character(len=80) mylnbuf boolinside=.false. boolhastag=.false. do while (.true.) read(ilun,'(80a)',end=1949) mylnbuf if (boolinside) then if (mylnbuf(1:14)=='!! BEGIN Amber') then write(iout,"('Error: Nested interface tags')") call mexit(iout,1) elseif (mylnbuf(1:14)=='!! END Amber') then boolinside=.false. else write(iout,'(80a)') mylnbuf(1:79) endif elseif (mylnbuf(1:14)=='!! BEGIN Amber') then write(iout,*) write(iout,"(' The Interface script used to generate the input file:')") write(iout,*) boolinside=.true. boolhastag=.true. elseif (mylnbuf(1:14)=='!! END Amber') then write(iout,"('Error: unclosed interface tag')") call mexit(iout,1) endif enddo 1949 rewind(ilun) if (boolhastag) write(iout,"(79('-'))") write(iout,*) write(iout,"(' Here is the input file:')") write(iout,*) boolinside=.false. do while (.true.) read(ilun,'(80a)',end=9562) mylnbuf if (boolinside) then if (mylnbuf(1:14)=='!! END Amber') then boolinside=.false. endif elseif (mylnbuf(1:14)=='!! BEGIN Amber') then boolinside=.true. else write(iout,'(80a)') mylnbuf(1:79) endif enddo 9562 rewind(ilun) return end subroutine myechoin subroutine mynmlsrc(srchkey,ilun,ifound) ! Author: Meng-Juei Hsieh implicit none character(len=*) srchkey integer ilun character(len=80) realkey character(len=80) mylnbuf integer ifound, ilen, ipos ilen = len(srchkey) realkey(2:ilen+1)=srchkey(1:ilen) ifound = 0 ipos = 0 do while (ifound == 0) read(ilun,'(80a)',end=824,err=9528) mylnbuf realkey(1:1)='&' ipos=index(mylnbuf, realkey(1:ilen+1), .false.) if (ipos>0) then ifound = 1 else realkey(1:1)='$' ipos=index(mylnbuf, realkey(1:ilen+1), .false.) if (ipos>0) ifound = 1 endif enddo 824 continue !write(6,*) realkey(1:ilen+1) if (ifound == 1 ) then backspace(ilun) else if (ifound == 0) then rewind(ilun) else write(6,*) "mynmlsrc Error: exception" call mexit(6,1) endif return 9528 rewind(ilun) return end subroutine mynmlsrc : TENTH = ONE/TEN double precision, parameter :: ELEVENTH = ONE/ELEVEN double precision, parameter :: TWELFTH = ONE/TWELVE ! dec cpp requires this comment line to avoid losing beginning spaces integer ilbopt, ilbnob, lbverb, lbfreq, dozero, plevel common/extra_int/ilbopt, ilbnob, lbverb, lbfreq, dozero, plevel double precision lbwght common/extra_real/lbwght logical master common/extra_logical/master !+ Specification and control of Amber's Input/Output ! File names character(len=512) groupbuffer character(len=80) mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe,redir(8), & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version character owrite common /files/ groupbuffer, mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe, owrite, & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version ! put this in a seperate common block to stop the compiler from ! complaining about misalignment integer numgroup common/nmgrp/ numgroup ! File units ! An I/O Unit resource manager does not exist. integer MDCRD_UNIT integer MDEN_UNIT integer MDINFO_UNIT integer MDVEL_UNIT parameter ( MDINFO_UNIT = 7 ) parameter ( MDCRD_UNIT = 12 ) parameter ( MDEN_UNIT = 15 ) parameter ( MDVEL_UNIT = 13 ) integer, parameter :: CNSTPH_UNIT = 18, CPOUT_UNIT = 19 ! 18 was picked because CNSTPH uses it; conflicts are not expected. !integer MMTSB_UNIT !parameter ( MMTSB_UNIT = 18 ) ! File related controls and options character(len=80) title,title1 common/runhed/ title, title1 logical mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version common/mdin_flags/mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version integer BC_HULP ! size in integers of common HULP parameter ( BC_HULP = 9 ) integer ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ,mdm_nstep & !#endif ; ! line terminator for free-form version common/hulp/ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ! this variable is last in the common and does not increment BC_HULP ! ,mdm_nstep & !#endif ; ! line terminator for free-form version ! NMRRDR : Contains information about input/output file redirection ! REDIR and IREDIR contain information regarding ! LISTIN, LISTOUT, READNMR, NOESY, SHIFTS, DUMPAVE, ! PCSHIFT and DIPOLE respectively. If IREDIR(I) > 0, ! then that input/output has been redirected. integer iredir(8) common/nmrrdr/redir,iredir ! 23320 is the number of reals in common RPARMS; 1200 is the ! number of ints in IPARMS. (Change these if you change the sizes below). double precision rk(5000),req(5000),tk(900),teq(900),pk(1200), & pn(1200),phase(1200),cn1(1830),cn2(1830),solty(60), & gamc(1200),gams(1200),fmn(1200), & asol(200),bsol(200),hbcut(200) common/rparms/rk,req,tk,teq,pk, & pn,phase,cn1,cn2,solty, & gamc,gams,fmn, & asol,bsol,hbcut integer ipn(1200) common/iparms/ipn ! NPHB is the number of h-bond parameters. NIMPRP is the number of ! improper torsional parameters (NPTRA-NIMPRP is the number of regular ! torsional parameters). integer numbnd,numang,nptra,nphb,nimprp common/prmlim/numbnd,numang,nptra,nphb,nimprp ! --- following should not need to be modified unless you are adding ! more variables to the "locmem" style of memory management ! 181 is the size of the MEMORY common block: ! MFC changing indices into IX and X to more rational names ! I12 = Iibh ! I14 = Ijbh ! I16 = Iicbh ! I18 = Iiba ! I20 = Ijba ! I22 = Iicba integer natom,nres,nbonh,nbona,ntheth,ntheta,nphih, & nphia,nnb,ntypes,nconp,maxmem,nwdvar,maxnb,nparm, & natc,ibelly,natbel,ishake,nmxrs, & mxsub,natyp,npdec,i02,i04,i06,i08,i10, & iibh,ijbh,iicbh,iiba,ijba,iicba, & i24,i26,i28,i30,i32,i34,i36,i38,i40, & i42,i44,i46,i48,i50,i52,i54,i56,i58,ibellygp, & icnstrgp,i64,i65,i66,i68, & i70,i72,i74,i76,i78,i79,i80,i82,i84,i86, & icpstinf,icpresst,icptrsct, icpptcnt, & l05,l10,l15,lwinv,lpol,lcrd,lforce,l36,lvel,lvel2,l45,l50, & lcrdr,l60,l65,lmass,l75,l80,l85,l90,l95,l96,l97,l98,lfrctmp, & l105,l110,l115,l120,l125,l130,l135,l140,l145,l150, & l165,l170,l175,l180,l185,l186,l187,l188,l189,l190, & lcpcrg,lcpene, & m02,m04,m06,m08,m10,m12,m14,m16,m18,i01, & lastr,lasti,lasth,lastpr,lastrst,lastist,ncopy ! lastpr,lastrst,lastist,nbper,ngper,ndper,ifpert,lpolp, ncopy, & ! imask1,imask2,numadjst,mxadjmsk ! 1 2 3 4 5 6 7 8 9 10 common/memory/ & natom ,nres ,nbonh ,nbona ,ntheth ,ntheta ,nphih , & ! 7 nphia ,nnb ,ntypes,nconp ,maxmem ,nwdvar ,maxnb ,nparm , & !15 natc ,ibelly,natbel,ishake,nmxrs , & !20 mxsub ,natyp ,npdec ,i02 ,i04 ,i06 ,i08 ,i10, & !28 iibh ,ijbh ,iicbh ,iiba ,ijba ,iicba , & !34 i24 ,i26 ,i28 ,i30 ,i32 ,i34 ,i36 ,i38 ,i40 , & !43 i42 ,i44 ,i46 ,i48 ,i50 ,i52 ,i54 ,i56 ,i58 ,ibellygp,& !53 icnstrgp,i64 ,i65 ,i66 ,i68 , & !58 i70 ,i72 ,i74 ,i76 ,i78 ,i79 ,i80 ,i82 , & !66 i84 ,i86 , & !68 icpstinf,icpresst,icptrsct, icpptcnt, & !72 l05 ,l10 ,l15 ,lwinv ,lpol ,lcrd ,lforce,l36 ,lvel ,lvel2 , & !82 l45 ,l50 , & !84 lcrdr ,l60 ,l65 ,lmass ,l75 ,l80 ,l85 ,l90 ,l95 ,l96 , & !94 l97 ,l98 ,lfrctmp, & !97 l105 ,l110 ,l115 ,l120 ,l125 ,l130 ,l135 ,l140 ,l145 ,l150 , & !107 l165 ,l170 ,l175 ,l180 ,l185 ,l186 ,l187 ,l188 ,l189 ,l190 , & !117 lcpcrg,lcpene, & !119 m02 ,m04 ,m06 ,m08 ,m10 ,m12 ,m14 ,m16 ,m18 ,i01 , & !129 lastr ,lasti ,lasth ,lastpr,lastrst,lastist,ncopy !136 !iifstwt,iifstwr,nrealb,nintb,nholb,npairb,lastr ,lasti ,lasth , & !138 !lastpr ,lastrst,lastist,nbper,ngper,ndper,ifpert,lpolp ,ncopy, & !147 !imask1 ,imask2 ,numadjst,mxadjmsk !151 ! common block sizes: integer bc_boxi,bc_boxr parameter(bc_boxi=7) parameter(bc_boxr=611) ! ... floats: double precision box,cut,dielc,rad,wel,radhb,welhb, & cutcap,xcap,ycap,zcap,fcap,rwell common/boxr/box(3),cut,dielc, & !5 cutcap,xcap,ycap,zcap,fcap,rwell, & !11 rad(100),wel(100),radhb(200),welhb(200) !611 ! ... integers: integer ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp common/boxi/ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp !-------------BEGIN md.h ------------------------------------------------ integer BC_MDI ! size in integers of common block mdi integer BC_MDR ! size in _REAL_s of common block mdr ! ... integer variables: integer nrp,nspm,ig,ntx,ntcx, &!5 ntxo,ntt,ntp,ntr,init, &!10 ntcm,nscm,isolvp,nsolut,klambda, &!15 ntc,ntcc,ntf,ntid,ntn, &!20 ntnb,nsnb,ndfmin,nstlim,nrc, &!25 ntrx,npscal,imin,maxcyc,ncyc, &!30 ntmin,irest,jfastw, &!33 ibgwat,ienwat,iorwat, &!36 iwatpr,nsolw,igb,iyammp, &!40 gbsa,vrand,iwrap,nrespa,irespa,nrespai,icfe, &!47 rbornstat,ivcap,iconstreff, &!50 idecomp,icnstph,ntcnstph,maxdup, &!54 ipb, inp, ibgion, ienion !58 parameter (BC_MDI=56) common/mdi/nrp,nspm,ig, & ntx,ntcx,ntxo,ntt,ntp,ntr,init,ntcm,nscm, & isolvp,nsolut,ntc,ntcc,ntf,ntid,ntn,ntnb,nsnb,ndfmin, & nstlim,nrc,ntrx,npscal,imin,maxcyc,ncyc,ntmin, & irest,jfastw,ibgwat,ienwat,iorwat, & iwatpr,nsolw,igb,iyammp,gbsa,vrand, & iwrap,nrespa,irespa,nrespai,icfe,rbornstat, & ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup, & ipb,inp,ibgion,ienion ! ... floats: double precision t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & !9 tol,taur,dx0,drms,vlimit,rbtarg(8),tmass,tmassinv, & !24 kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & !31 cut_inner,clambda,saltcon, & !34 solvph,rgbmax,fsmax,restraint_wt !38 parameter (BC_MDR=38) common/mdr/t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & tol,taur,dx0,drms,vlimit,rbtarg,tmass,tmassinv, & kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & cut_inner,clambda,saltcon, & solvph,rgbmax,fsmax,restraint_wt ! ... strings: character(len=4) iwtnm,iowtnm,ihwtnm character(len=256) restraintmask,bellymask common/mds/ restraintmask,bellymask,iwtnm,iowtnm,ihwtnm(2) !-------------END md.h ------------------------------------------------ close(unit=8) ! ----- SET THE DEFAULT VALUES FOR SOME VARIABLES ----- nrp = natom if (ifbox == 1) write(6, '(/5x,''BOX TYPE: RECTILINEAR'')') if (ifbox == 2) write(6, '(/5x,''BOX TYPE: TRUNCATED OCTAHEDRON'')') if (ifbox == 3) write(6, '(/5x,''BOX TYPE: GENERAL'')') nsolut = nrp if ( nscm > 0 .and. ntb == 0 ) then ndfmin = 6 ! both translation and rotation com motion removed if (nsolut == 1) ndfmin = 3 if (nsolut == 2) ndfmin = 5 else if ( nscm > 0 ) then ndfmin = 3 ! just translation com will be removed else ndfmin = 0 end if if (ibelly > 0) then ! No COM Motion Removal, ever. nscm = 9999999 ndfmin = 0 end if if(nscm <= 0) nscm = 9999999 init = 3 if (irest > 0) then ! init = 4 write(6,*) "We are sorry, but irest have to be 0" call mexit(6,0) endif !if (scnb == ZERO ) scnb = TWO if (dielc <= ZERO ) dielc = ONE if (tautp <= ZERO ) tautp = 0.2d0 if (taup <= ZERO ) taup = 0.2d0 ! RESET THE CAP IF NEEDED if(ivcap == 2) ifcap = 0 ! PRINT DATA CHARACTERIZING THE RUN nmropt = 0 iesp = 0 ipol = 0 write(6,9328) write(6,9008) title write(6,'(/a)') 'General flags:' write(6,'(5x,2(a,i8))') 'imin =',imin,', nmropt =',nmropt write(6,'(/a)') 'Nature and format of input:' write(6,'(5x,4(a,i8))') 'ntx =',ntx,', irest =',irest, & ', ntrx =',ntrx write(6,'(/a)') 'Nature and format of output:' write(6,'(5x,4(a,i8))') 'ntxo =',ntxo,', ntpr =',ntpr, & ', ntrx =',ntrx,', ntwr =',ntwr write(6,'(5x,4(a,i8))') 'iwrap =',iwrap,', ntwx =',ntwx, & ', ntwv =',ntwv,', ntwe =',ntwe write(6,'(5x,3(a,i8),a,i7)') 'ioutfm =',ioutfm, & ', ntwprt =',ntwprt, & ', idecomp =',idecomp,', rbornstat=',rbornstat write(6,'(/a)') 'Potential function:' write(6,'(5x,5(a,i8))') 'ntf =',ntf,', ntb =',ntb, & ', igb =',igb,', nsnb =',nsnb write(6,'(5x,4(a,i8))') 'ipol =',ipol,', gbsa =',gbsa, & ', iesp =',iesp write(6,'(5x,3(a,f10.5))') 'dielc =',dielc, & ', cut =',cut,', intdiel =',intdiel !write(6,'(5x,3(a,f10.5))') 'scnb =',scnb, & ! ', scee =',scee write(6,'(/a)') 'Frozen or restrained atoms:' write(6,'(5x,4(a,i8))') 'ibelly =',ibelly,', ntr =',ntr if( imin /= 0 ) then write(6,'(/a)') 'Energy minimization:' ! print inputable variables applicable to all minimization methods. write(6,'(5x,4(a,i8))') 'maxcyc =',maxcyc,', ncyc =',ncyc, & ', ntmin =',ntmin write(6,'(5x,2(a,f10.5))') 'dx0 =',dx0, ', drms =',drms ! ! Input flag ntmin determines the method of minimization ! select case ( ntmin ) ! case ( 0, 1, 2 ) ! ! no specific output ! case default ! ! invalid ntmin ! write(6,'(a,i3)') ' ERROR: Invalid NTMIN.', ntmin ! call mexit(6, 1) ! end select else write(6,'(/a)') 'Molecular dynamics:' write(6,'(5x,4(a,i8))') 'nstlim =',nstlim,', nscm =',nscm, & ', nrespa =',nrespa write(6,'(5x,3(a,f10.5))') 't =',t, & ', dt =',dt,', vlimit =',vlimit end if cut = cut*cut cut_inner = cut_inner*cut_inner ! If user has requested Poisson-Boltzmann electrostatics, set up variables call pb_init(ifcap,natom,nres,ntypes,nbonh,nbona,ix(i02),ix(i04),ix(i06),ix(i08),ix(i10),& ix(iibh),ix(ijbh),ix(iiba),ix(ijba),ix(ibellygp),ih(m02),ih(m04),ih(m06),x(l15),x(l97)) ntmp = 0; ngrp = 0 if(idecomp > 0) then write(6,9428) call rgroup(natom,ntmp,nres,ngrp,ix(i02),ih(m02), & ih(m04),ih(m06),ih(m08),ix(ibellygp), & jgroup,index,irespw,npdec, & x(l60),x(lcrdr),.false.,.false.,.false.,idecomp,5,.true.) end if ! checking of not supported options if( ifcap /= 0 ) then write(6, '(a,i3)') ' ERROR: cap water not supported.', ifcap call mexit(6, 1) end if if( igb /= 10) then write(6, '(a,i3)') ' ERROR: Non PB potential not supported', igb call mexit(6, 1) end if ! --- checks on bogus data --- inerr = 0 if (imin < 0 .or. imin > 1) then write(6,'(/2x,a,i3,a)') 'IMIN (',imin,') must be 0 or 1.' inerr = 1 end if if (ntx < 1 .or. ntx > 7) then write(6,'(/2x,a,i3,a)') 'NTX (',ntx,') must be in 1..7' inerr = 1 end if ! ---WARNINGS: if (inerr == 1) then write(6, '(/,a)') ' *** input error(s)' call mexit(6,1) end if ! DEBUG input; force checking call load_debug(5) ! Standard format statements: 9328 format(/80('-')/,' 2. CONTROL DATA FOR THE RUN',/80('-')/) 9428 format(/4x,'LOADING THE DECOMP ATOMS AS GROUPS',/) 9008 format(a80) end subroutine mdread2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Emit defined preprocessor names, ie, flags. subroutine printflags() implicit none integer max_line_length parameter ( max_line_length = 80 ) character(len=max_line_length) line ! output string of active flags integer n ! len(line) line = '| Flags:' n = 8 call printflags2(' ',1,n,line,.true.) return end subroutine printflags !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Primitive pre-Fortran90 implementation of printflags. subroutine printflags2(flag,flag_len,line_len,line,last) implicit none integer max_line_length parameter ( max_line_length = 80 ) character(*) flag ! flag name with blank prefix, intent(in) integer flag_len ! len(flag), intent(in) integer line_len ! len(line), intent(inout) character(len=max_line_length) line ! intent(inout) logical last ! is this the last flag ?, intent(in) if (line_len + flag_len > max_line_length) then write( 6,'(a)') line ! begin another line line = '| Flags:' line_len=8 end if line=line(1:line_len) // flag(1:flag_len) line_len=line_len+flag_len if(last)write( 6,'(a)') line return end subroutine printflags2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Check the range of a float; abort on illegal values. subroutine float_legal_range(string,param,lo,hi) implicit none double precision param,lo,hi character(len=*)string if ( param < lo .or. param > hi )then write(6,59) write(6,60)string,param write(6,61) write(6,62)lo,hi write(6,63) call mexit(6,1) end if 59 format(/,1x,'Ewald PARAMETER RANGE CHECKING: ') 60 format(1x,'parameter ',a,' has value ',e12.5) 61 format(1x,'This is outside the legal range') 62 format(1x,'Lower limit: ',e12.5,' Upper limit: ',e12.5) 63 format(1x,'Check ew_legal.h') return end subroutine float_legal_range !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Check the range of an integer; abort on illegal values. subroutine int_legal_range(string,param,lo,hi) implicit none integer param,lo,hi character(len=*)string if ( param < lo .or. param > hi )then write(6,59) write(6,60)string,param write(6,61) write(6,62)lo,hi write(6,63) call mexit(6,1) end if 59 format(/,1x,'PARAMETER RANGE CHECKING: ') 60 format(1x,'parameter ',a,' has value ',i8) 61 format(1x,'This is outside the legal range') 62 format(1x,'Lower limit: ',i8,' Upper limit: ',i8) 63 format(1x,'The limits may be adjustable; search in the .h files ') return end subroutine int_legal_range !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Check the range of an integer option; abort on illegal values. subroutine opt_legal_range(string,param,lo,hi) implicit none integer param,lo,hi character(len=*)string if ( param < lo .or. param > hi )then write(6,59) write(6,60)string,param write(6,61) write(6,62)lo,hi write(6,63) call mexit(6,1) end if 59 format(/,1x,'Ewald OPTION CHECKING: ') 60 format(1x,'option ',a,' has value ',i5) 61 format(1x,'This is outside the legal range') 62 format(1x,'Lower limit: ',i5,' Upper limit: ',i5) 63 format(1x,'Check the manual') return end subroutine opt_legal_range !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ PBMD module parsing and initialization subroutine pb_read ! Module variables use poisson_boltzmann use dispersion_cavity use solvent_accessibility implicit none ! Common variables !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+pb_def.h 1 cpp macro-definitions ! maximum number of neighboring atoms per atom ! maximum levels of focusing in fdpb ! maximum number of grid charges for boundary conditions !+End of pb_def.h !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Specification and control of Amber's Input/Output ! File names character(len=512) groupbuffer character(len=80) mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe,redir(8), & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version character owrite common /files/ groupbuffer, mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmr, mincor, & vecs, radii, freqe, owrite, & rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt & !#ifdef MMTSB ! ,mmtsb_setup_file & !#endif !cnt N.TAKADA: For MDM !#ifdef MDM_PDB ! ,mdpdb & !#endif !cnt N.TAKADA: For MDM ; ! line terminator for free-form version ! put this in a seperate common block to stop the compiler from ! complaining about misalignment integer numgroup common/nmgrp/ numgroup ! File units ! An I/O Unit resource manager does not exist. integer MDCRD_UNIT integer MDEN_UNIT integer MDINFO_UNIT integer MDVEL_UNIT parameter ( MDINFO_UNIT = 7 ) parameter ( MDCRD_UNIT = 12 ) parameter ( MDEN_UNIT = 15 ) parameter ( MDVEL_UNIT = 13 ) integer, parameter :: CNSTPH_UNIT = 18, CPOUT_UNIT = 19 ! 18 was picked because CNSTPH uses it; conflicts are not expected. !integer MMTSB_UNIT !parameter ( MMTSB_UNIT = 18 ) ! File related controls and options character(len=80) title,title1 common/runhed/ title, title1 logical mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version common/mdin_flags/mdin_pb & !#ifdef MDM_MD ! ,mdin_mdm & !#endif !#ifdef MDM_PDB ! ,mdin_pdb & !#endif ; ! line terminator for free-form version integer BC_HULP ! size in integers of common HULP parameter ( BC_HULP = 9 ) integer ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ,mdm_nstep & !#endif ; ! line terminator for free-form version common/hulp/ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave & !#ifdef MDM_MD ! ! this variable is last in the common and does not increment BC_HULP ! ,mdm_nstep & !#endif ; ! line terminator for free-form version ! NMRRDR : Contains information about input/output file redirection ! REDIR and IREDIR contain information regarding ! LISTIN, LISTOUT, READNMR, NOESY, SHIFTS, DUMPAVE, ! PCSHIFT and DIPOLE respectively. If IREDIR(I) > 0, ! then that input/output has been redirected. integer iredir(8) common/nmrrdr/redir,iredir !-------------BEGIN md.h ------------------------------------------------ integer BC_MDI ! size in integers of common block mdi integer BC_MDR ! size in _REAL_s of common block mdr ! ... integer variables: integer nrp,nspm,ig,ntx,ntcx, &!5 ntxo,ntt,ntp,ntr,init, &!10 ntcm,nscm,isolvp,nsolut,klambda, &!15 ntc,ntcc,ntf,ntid,ntn, &!20 ntnb,nsnb,ndfmin,nstlim,nrc, &!25 ntrx,npscal,imin,maxcyc,ncyc, &!30 ntmin,irest,jfastw, &!33 ibgwat,ienwat,iorwat, &!36 iwatpr,nsolw,igb,iyammp, &!40 gbsa,vrand,iwrap,nrespa,irespa,nrespai,icfe, &!47 rbornstat,ivcap,iconstreff, &!50 idecomp,icnstph,ntcnstph,maxdup, &!54 ipb, inp, ibgion, ienion !58 parameter (BC_MDI=56) common/mdi/nrp,nspm,ig, & ntx,ntcx,ntxo,ntt,ntp,ntr,init,ntcm,nscm, & isolvp,nsolut,ntc,ntcc,ntf,ntid,ntn,ntnb,nsnb,ndfmin, & nstlim,nrc,ntrx,npscal,imin,maxcyc,ncyc,ntmin, & irest,jfastw,ibgwat,ienwat,iorwat, & iwatpr,nsolw,igb,iyammp,gbsa,vrand, & iwrap,nrespa,irespa,nrespai,icfe,rbornstat, & ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup, & ipb,inp,ibgion,ienion ! ... floats: double precision t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & !9 tol,taur,dx0,drms,vlimit,rbtarg(8),tmass,tmassinv, & !24 kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & !31 cut_inner,clambda,saltcon, & !34 solvph,rgbmax,fsmax,restraint_wt !38 parameter (BC_MDR=38) common/mdr/t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & tol,taur,dx0,drms,vlimit,rbtarg,tmass,tmassinv, & kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & cut_inner,clambda,saltcon, & solvph,rgbmax,fsmax,restraint_wt ! ... strings: character(len=4) iwtnm,iowtnm,ihwtnm character(len=256) restraintmask,bellymask common/mds/ restraintmask,bellymask,iwtnm,iowtnm,ihwtnm(2) !-------------END md.h ------------------------------------------------ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ pb_md.h PB variables shared with calling routines logical pbverbose logical pbprint logical pbgrid logical pbinit common /pb_mdl/ pbverbose common /pb_mdl/ pbprint common /pb_mdl/ pbgrid common /pb_mdl/ pbinit integer npbstep integer nsaslag integer npbgrid integer nsnbr integer nsnba integer ntnbr integer ntnba integer dofd integer ndofd integer dosas integer ndosas common /pb_mdi/ npbstep common /pb_mdi/ nsaslag common /pb_mdi/ npbgrid common /pb_mdi/ nsnbr common /pb_mdi/ nsnba common /pb_mdi/ ntnbr common /pb_mdi/ ntnba common /pb_mdi/ dofd common /pb_mdi/ ndofd common /pb_mdi/ dosas common /pb_mdi/ ndosas !+ End of pb_md.h !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Local variables integer npbverb, l, phiout, scalec double precision space ! begin code namelist /pb/ epsin, epsout, smoothopt, istrng, pbtemp, & radiopt, dprob, iprob, npbopt, solvopt, accept, maxitn, & fillratio, space, nbuffer, nfocus, fscale, npbgrid, & arcres,dbfopt,bcopt,scalec,eneopt,frcopt,cutres,cutfd, & cutnb, nsnbr, nsnba,phiout, phiform, npbverb, npopt, & decompopt, use_rmin, sprob, vprob, rhow_effect, use_sav, & cavity_surften, cavity_offset, maxsph, maxarc, & cutsa, ndofd, ndosas, fmiccg, ivalence, laccept, wsor, & lwsor, pbkappa, radinc, expthresh, offx, offy, offz, & sepbuf, mpopt, lmax, maxarcdot ! initialize with default values outphi = .false. phiout = 0 phiform = 0 epsin = 1.0d0 epsout = 80.0d0 istrng = 0.0d0 ivalence = 1.0d0 pbtemp = 300.0d0 scalec = 0 scalerf = .false. srsas = .true. smoothopt = 0 radiopt = 1 npopt = 2 decompopt = 1 use_rmin = 0 use_sav = 0 rhow_effect = 1.0d0 maxsph = 400 maxarcdot= 0 maxarc = 256 nbuffer = 0 dprob = 0.00d0 sprob = 1.60d0 vprob = 1.28d0 iprob = 2.00d0 radinc = sprob*0.5d0 expthresh = 0.2d0 arcres = 0.0625d0 cavity_surften = 0.04356d0 cavity_offset = -1.008d0 nfocus = 2 fscale = 8 level = 1 bcopt = 5 space = 0.5d0 savxm(1) = 0 savym(1) = 0 savzm(1) = 0 savxmym(1) = 0 savxmymzm(1) = 0 savh(1) = 0 savbcopt(1) = 0 savxm(nfocus) = 0 savym(nfocus) = 0 savzm(nfocus) = 0 savxmym(nfocus) = 0 savxmymzm(nfocus) = 0 savh(nfocus) = 0 savbcopt(nfocus) = 0 offx = 0.0d0 offy = 0.0d0 offz = 0.0d0 fillratio= 2.0d0 npbopt = 0 solvopt = 1 fmiccg = -0.30d0 accept = 1.0d-3 laccept = 0.1d0 wsor = 1.9d0 lwsor = 1.95d0 maxitn = 100 pbverbose = .false. pbprint = .true. pbgrid = .true. pbinit = .true. npbverb = 0 npbgrid = 1 ndofd = 1 dofd = 1 donpsa = .true. ndosas = 1 nsaslag = 100 nsnbr = 1 nsnba = 1 ntnbr = 1 ntnba = 1 cutres = 99.0d0 cutfd = 5.0d0 cutnb = 0.0d0 cutsa = 9.0d0 lastp = 0 pbgamma_int = 1.0 pbgamma_ext = 65.0 sepbuf = 4.0d0 mpopt = 0 lmax = 80 dbfopt = -1 eneopt = -1 frcopt = 0 ! reading parameters if ( mdin_pb ) then rewind 5 read(5, nml=pb) end if ! check and post processing of options if ( npbverb > 0 ) pbverbose = .true. dofd = ndofd dosas = ndosas epsin = epsin*eps0 epsout = epsout*eps0 istrng = fioni * istrng pbkappa = SQRT( 2.0d0 * istrng / (epsout * pbkb * pbtemp) ) ! set up solvent probes for different solvation components ! if end user requests different probes for different solvation ! components, no need to set up different probes here. if ( dprob == 0.0d0 ) dprob = sprob ! set up solvent accessible arc resolutions and limits if ( maxarcdot == 0) maxarcdot = 1200*nint(0.5d0/arcres) ! check dprob and iprob if ( dprob >= iprob ) then write(6, *) 'PB Bomb in pb_read(): solvent probe cannot be smaller than ion prob' call mexit(6, 1) end if if ( dprob <= space .and. ipb == 1 ) then write(6, *) 'PB Warning in pb_read(): setting grid spacing larger than solvent probe' write(6, *) 'may cause numerical instability if ipb=1' call mexit(6, 1) end if if ( dprob <= space .and. ipb == 2 ) then write(6, *) 'PB Bomb in pb_read(): solvent probe cannot be smaller than grid spacing if ipb=2' call mexit(6, 1) end if ! check pb options if ( igb == 10 .and. ipb /= 1 ) then write(6, *) 'PB Info in pb_read(): igb has been overwritten by ipb' endif ! check nonpolar options if ( inp /= npopt ) then write(6, *) 'PB Info in pb_read(): npopt has been overwritten with inp' npopt = inp endif if ( inp == 2 .and. radiopt == 0 ) then 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) end if if ( radiopt == 0 ) then donpsa = .false. end if ! check force options if ( eneopt == -1 ) then if ( dbfopt == 0 ) then eneopt = 1 else if ( dbfopt == 1 .or. dbfopt == -1 ) then eneopt = 2 else write(6, *) 'PB Bomb in pb_read(): only dbfopt= 0 or 1 are supported. dbfopt is replaced by eneopt' call mexit(6, 1) end if else if ( dbfopt == 0 .and. eneopt == 1 ) then write(6, *) 'PB Info in pb_read(): dbfopt has been overwritten by eneopt' else if ( dbfopt == 1 .and. eneopt == 2 ) then write(6, *) 'PB Info in pb_read(): dbfopt has been overwritten by eneopt' else if ( dbfopt /= -1 ) then write(6, *) 'PB Bomb in pb_read(): dbfopt is ignored' end if end if if ( eneopt < 1 .or. eneopt > 2 ) then write(6, *) 'PB Bomb in pb_read(): only eneopt= 1 or 2 are supported' call mexit(6, 1) end if if ( frcopt < 0 .or. eneopt > 3 ) then write(6, *) 'PB Bomb in pb_read(): only frcopt= 0-3 are supported' call mexit(6, 1) end if if ( eneopt == 1 .and. ( frcopt == 2 .or. frcopt == 3 ) ) then write(6, *) 'PB Bomb in pb_read(): this combination of eneopt and frcopt is not supported' call mexit(6, 1) end if if ( eneopt == 2 .and. frcopt == 1 ) then write(6, *) 'PB Bomb in pb_read(): this combination of eneopt and frcopt is not supported' call mexit(6, 1) end if if ( ( frcopt == 2 .or. frcopt == 3 ) .and. bcopt == 5 ) then bcopt = 6 write(6, *) 'PB Info in pb_read(): bcopt has been reset from 5 to 6 for frcopt= 2 or 3' end if if ( frcopt == 1 .and. bcopt == 6 ) then bcopt = 5 write(6, *) 'PB Info in pb_read(): bcopt has been reset from 6 to 5 for frcopt=1' end if if ( frcopt > 0 .and. smoothopt == 0 ) then smoothopt = 1 write(6, *) 'PB Info in pb_read(): smoothopt has been reset to 1 for force compuation' end if if ( ivcap /= 0 .and. (eneopt /= 1 .or. frcopt /= 0) ) then write(6,*) "cap water should only be run with eneopt=1 and frcopt=0" call mexit(6,1) endif ! check numerical solution options if ( npbopt == 0 .and. solvopt > 4 ) then write(6, *) 'PB Bomb in pb_read(): solvopt>4 cannot be used to solve linear PB equation' call mexit(6, 1) endif if ( solvopt == 2 .and. bcopt == 6 ) then write(6, *) 'PB Bomb in pb_read(): bcopt=6 cannot be used with solvopt=2' call mexit(6, 1) end if if ( npbopt == 1 .and. solvopt > 6 ) then write(6, *) 'PB Bomb in pb_read(): nonsupported solvopt (>6) for e nonlinear PB equation' call mexit(6, 1) endif if ( npbopt == 1 .and. eneopt == 2 ) then eneopt = 1 write(6, *) 'PB Info in pb_read(): eneopt has been reset to be 1 for nonlinear PB equation' end if if ( npbopt == 1 .and. frcopt > 0 ) then write(6, *) 'PB Bomb in pb_read(): force computatoin is not supported for nonlinear PB equation' call mexit(6, 1) end if ! check cutoff options if ( cutfd > cutsa ) then cutsa = cutfd end if if ( cutnb /= 0 .and. cutfd > cutnb ) then cutnb = cutfd write(6, *) 'PB Info in pb_read(): cutnb has been reset to be equal to cutfd' end if if ( cutnb /= 0 .and. bcopt == 6 ) then cutnb = 0 write(6, *) 'PB Info in pb_read(): cutnb has been reset to 0 with bcopt=6', cutnb, bcopt end if if ( max(cutnb,cutsa,cutfd) > cutres ) then write(6, *) 'PB Bomb in pb_read(): cutnb/cutfd must be <= cutres', cutnb, cutfd, cutres call mexit(6, 1) end if if ( cutnb == 0 .and. eneopt == 1 .and. bcopt /= 6 ) then write(6, *) 'PB Bomb in pb_read(): cutnb=0 cannot be used with eneopt=1', cutnb, eneopt call mexit(6, 1) end if cutfd = cutfd**2 cutsa = cutsa**2 cutnb = cutnb**2 cutres = cutres**2 ! set buffer zone between the fine FD grid boundary and the solute surface: if ( nbuffer == 0 ) then if ( istrng == 0.0d0 ) then nbuffer = int(2.0d0*sprob/space)+1 else if ( sprob >= iprob ) then nbuffer = int(2.0d0*sprob/space)+1 else nbuffer = int(2.0d0*iprob/space)+1 end if end if if ( nbuffer >= fscale ) then nbuffer = 2*nbuffer+1 else nbuffer = 2*fscale+1 end if end if ! set flag to scale induced surface charges: if ( scalec == 1) scalerf = .true. ! set saved grid options savbcopt(1) = bcopt do l = 2, nfocus savbcopt(l) = 0 end do savh(nfocus) = space do l = nfocus - 1, 1, -1 savh(l) = savh(l+1)*fscale end do do l = 1, nfocus savxm(l) = 1 savym(l) = 1 savzm(l) = 1 savxmym(l) = 1 savxmymzm(l) = 1 end do ! set phimap output options when requested if ( phiout == 1 ) then outphi = .true. radiopt = 2 donpsa = .false. npopt = 0 end if end subroutine pb_read subroutine myechoin (ilun,iout) ! Rewritten by: Meng-Juei Hsieh ! A rewrite of echoin implicit none integer ilun, iout logical boolinside, boolhastag character(len=80) mylnbuf boolinside=.false. boolhastag=.false. do while (.true.) read(ilun,'(80a)',end=1949) mylnbuf if (boolinside) then if (mylnbuf(1:14)=='!! BEGIN Amber') then write(iout,"('Error: Nested interface tags')") call mexit(iout,1) elseif (mylnbuf(1:14)=='!! END Amber') then boolinside=.false. else write(iout,'(80a)') mylnbuf(1:79) endif elseif (mylnbuf(1:14)=='!! BEGIN Amber') then write(iout,*) write(iout,"(' The Interface script used to generate the input file:')") write(iout,*) boolinside=.true. boolhastag=.true. elseif (mylnbuf(1:14)=='!! END Amber') then write(iout,"('Error: unclosed interface tag')") call mexit(iout,1) endif enddo 1949 rewind(ilun) if (boolhastag) write(iout,"(79('-'))") write(iout,*) write(iout,"(' Here is the input file:')") write(iout,*) boolinside=.false. do while (.true.) read(ilun,'(80a)',end=9562) mylnbuf if (boolinside) then if (mylnbuf(1:14)=='!! END Amber') then boolinside=.false. endif elseif (mylnbuf(1:14)=='!! BEGIN Amber') then boolinside=.true. else write(iout,'(80a)') mylnbuf(1:79) endif enddo 9562 rewind(ilun) return end subroutine myechoin subroutine mynmlsrc(srchkey,ilun,ifound) ! Author: Meng-Juei Hsieh implicit none character(len=*) srchkey integer ilun character(len=80) realkey character(len=80) mylnbuf integer ifound, ilen, ipos ilen = len(srchkey) realkey(2:ilen+1)=srchkey(1:ilen) ifound = 0 ipos = 0 do while (ifound == 0) read(ilun,'(80a)',end=824,err=9528) mylnbuf realkey(1:1)='&' ipos=index(mylnbuf, realkey(1:ilen+1), .false.) if (ipos>0) then ifound = 1 else realkey(1:1)='$' ipos=index(mylnbuf, realkey(1:ilen+1), .false.) if (ipos>0) ifound = 1 endif enddo 824 continue !write(6,*) realkey(1:ilen+1) if (ifound == 1 ) then backspace(ilun) else if (ifound == 0) then rewind(ilun) else write(6,*) "mynmlsrc Error: exception" call mexit(6,1) endif return 9528 rewind(ilun) return end subroutine mynmlsrc