!************************************************************************ ! AMBER ** ! ** ! Copyright (c) 2003 ** ! Regents of the University of California ** ! All Rights Reserved. ** ! ** ! This software provided pursuant to a license agreement containing ** ! restrictions on its disclosure, duplication, and use. This software ** ! contains confidential and proprietary information, and may not be ** ! extracted or distributed, in whole or in part, for any purpose ** ! whatsoever, without the express written permission of the authors. ** ! This notice, and the associated author list, must be attached to ** ! all copies, or extracts, of this software. Any additional ** ! restrictions set forth in the license agreement also apply to this ** ! software. ** !************************************************************************ !+ Specification and control of Amber's working precision ! Description: ! Preprocessor directives that characterize the floating-point ! working precision as single or double precision. ! The current scheme guarantees internal consistency at the expense ! of flexibility. A need for flexibility has yet to appear. ! The preprocessor guard prevents multiple, and thus ! inconsistent, definitions. ! The default working precision is double precision. ! User control of the working precision at build time should be ! exercised via the preprocessor name _REAL_. ! To build a single precision Amber use ! make -e AMBERBUILDFLAGS=' -D_REAL_ ' ! The preprocessor names that characterize the precision are ! _REAL_ precision type specifier. ! Use _REAL_ foo as the precision independent ! notation for double precision foo and real foo. ! AMBER_MPI_REAL ! MPI precision type specifier. ! Use AMBER_MPI_REAL as the precision independent ! notation for MPI_DOUBLE_PRECISION and MPI_REAL. ! D_OR_S() precision prefix for the BLAS and LAPACK Library routines. ! Use, e.g., D_OR_S()axpy(...) as the precision independent ! notation for daxpy(...) and saxpy(...). ! DPREC defined when the working precision is double; ! undefined when the working precision is single. ! VD_OR_VS() precision prefix for the Intel Vector Math Library routines. ! Use, e.g., VD_OR_VS()exp(...) as the precision independent ! notation for vdexp(...) and vsexp(...). !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine rdparm1 here] subroutine rdparm1(nf) use parms, only: numbnd,numang,nptra,nphb,MAX_ATOM_TYPE,MAX_BOND_TYPE, & charmm use molecule implicit none !+ Specification and control of Amber's working precision !-------------BEGIN md.h ------------------------------------------------ integer BC_MDI ! size in integers of common block mdi integer BC_MDR ! size in Reals 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,alpb,iyammp, &!41 gbsa,vrand,iwrap,nrespa,irespa,nrespai,icfe, &!48 rbornstat,ivcap,iconstreff, &!51 neb,vv,tmode,ipol,iesp,ievb,nodeid,num_noshake, &!59 idecomp,icnstph,ntcnstph,maxdup,numexchg,repcrd,numwatkeep,hybridgb, &!67 ibgion,ienion,profile_mpi !70 parameter (BC_MDI=70) 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,alpb,iyammp,gbsa,vrand,numexchg,repcrd,numwatkeep,hybridgb, & iwrap,nrespa,irespa,nrespai,icfe,rbornstat, & ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup,neb,vv, & tmode,ipol,iesp,ievb,nodeid,num_noshake,ibgion,ienion, profile_mpi ! ... floats: double precision t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & !9 tol,taur,dx0,drms,vlimit,rbtarg(9),tmass,tmassinv, & !25 kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & !32 gbalpha,gbbeta,gbgamma,cut_inner,clambda,saltcon, & !38 solvph,rgbmax,fsmax,restraint_wt, & !42 skmin,skmax,vfac,gbneckscale,v11,v12,v22,kevb,evbt,Arad !52 parameter (BC_MDR=52) 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, & gbalpha,gbbeta,gbgamma,cut_inner,clambda,saltcon, & solvph,rgbmax,fsmax,restraint_wt,skmin,skmax,vfac,gbneckscale, & v11,v12,v22,kevb,evbt,Arad ! ... strings: character(len=4) iwtnm,iowtnm,ihwtnm character(len=256) restraintmask,bellymask,tgtfitmask,tgtrmsmask,noshakemask,crgmask common/mds/ restraintmask,bellymask,tgtfitmask,tgtrmsmask,noshakemask,crgmask, & iwtnm,iowtnm,ihwtnm(2) !-------------END md.h ------------------------------------------------ ! --- following should not need to be modified unless you are adding ! more variables to the "locmem" style of memory management ! 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,nparm, & natc,nattgtfit,nattgtrms,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,itgtfitgp,itgtrmsgp,i64,i65,i68, & i70,i72,i74,i76,i78,i79,i80,i82,i84,i86, & icpstinf,icpresst,icptrsct, icpptcnt, & l15,lwinv,lpol,lcrd,lforce,l36,lvel,lvel2,l45,l50, & lcrdr,l60,l65,lmass,l75,l80,l85,l90,l95,l96,l97,l98,l99,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, & iifstwt,iifstwr,nrealb,nintb,nholb,npairb,lastr,lasti,lasth, & lastpr,nbper,ngper,ndper,ifpert,ncopy, & imask1,imask2,numadjst,mxadjmsk,icphidx,icptpair,lfsg,lvsg,noshake ! 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,nparm , & !14 natc ,nattgtfit,nattgtrms,ibelly ,natbel,ishake,nmxrs , & !21 mxsub ,natyp ,npdec ,i02 ,i04 ,i06 ,i08 ,i10 , & !29 iibh ,ijbh ,iicbh ,iiba ,ijba ,iicba , & !35 i24 ,i26 ,i28 ,i30 ,i32 ,i34 ,i36 ,i38 ,i40 , & !44 i42 ,i44 ,i46 ,i48 ,i50 ,i52 ,i54 ,i56 ,i58 ,ibellygp,& !54 icnstrgp,itgtfitgp,itgtrmsgp,i64 ,i65 ,i68 , & !60 i70 ,i72 ,i74 ,i76 ,i78 ,i79 ,i80 ,i82 , & !68 i84 ,i86 , & !70 icpstinf,icpresst ,icptrsct ,icpptcnt, & !74 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 ,l99 ,lfrctmp , & !98 l105 ,l110 ,l115 ,l120 ,l125 ,l130 ,l135 ,l140 ,l145 ,l150 , & !108 l165 ,l170 ,l175 ,l180 ,l185 ,l186 ,l187 ,l188 ,l189 ,l190 , & !118 lcpcrg ,lcpene , & !120 m02 ,m04 ,m06 ,m08 ,m10 ,m12 ,m14 ,m16 ,m18 ,i01 , & !130 iifstwt ,iifstwr ,nrealb ,nintb ,nholb ,npairb,lastr ,lasti ,lasth , & !139 lastpr ,nbper ,ngper ,ndper ,ifpert, ncopy , & !145 imask1 ,imask2 ,numadjst ,mxadjmsk,icphidx,icptpair,lfsg,lvsg,noshake !154 ! BC_MEMORY is the size of the MEMORY common block !+ Specification and control of Amber's Input/Output ! File names character(len=4096) groupbuffer character(len=256) mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmrf, mincor, & vecs, radii, freqe,redir(9),rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt, evbin, evbout, mmtsb_setup_file,pimdout, & inptraj character owrite, facc common /files/ groupbuffer, mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmrf, mincor, & vecs, radii, freqe, owrite, facc,rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt, evbin, evbout, mmtsb_setup_file,pimdout, & inptraj ! put this in a separate common block to stop the compiler from ! complaining about misalignment integer numgroup, nslice common/nmgrp/ numgroup, nslice ! File units ! An I/O Unit resource manager does not exist. integer MDCRD_UNIT integer INPTRAJ_UNIT integer MDEN_UNIT integer MDINFO_UNIT integer MDVEL_UNIT parameter ( MDINFO_UNIT = 7 ) parameter ( MDCRD_UNIT = 12 ) parameter ( INPTRAJ_UNIT = 24 ) 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 ) !! !! EVB I/O unit !! integer, parameter :: evb_unit = 75 integer, parameter :: schlegel_unit = 80 !! FULL PIMD I/O unit integer, parameter :: pimd_unit = 277 ! File related controls and options character(len=80) title,title1 common/runhed/ title, title1 logical mdin_ewald,mdin_pb,mdin_amoeba common/mdin_flags/mdin_ewald,mdin_pb,mdin_amoeba integer BC_HULP ! size in integers of common HULP parameter ( BC_HULP = 9 ) integer ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave common/hulp/ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave ! 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(9) common/nmrrdr/redir,iredir ! common block sizes: ! ... floats: !+ Specification and control of Amber's working precision double precision box,cut,scnb,scee,dielc,rad,wel,radhb,welhb, & cutcap,xcap,ycap,zcap,fcap, & xlorth,ylorth,zlorth,xorth,yorth,zorth,forth, & rwell,xbox0,ybox0,zbox0 common/boxr/box(3),cut,scnb,scee,dielc,xbox0,ybox0,zbox0, & cutcap,xcap,ycap,zcap,fcap, & xlorth,ylorth,zlorth,xorth,yorth,zorth,forth, & rwell,rad(100),wel(100),radhb(200),welhb(200) ! ... integers: integer ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp common/boxi/ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp double precision extraboxdim parameter (extraboxdim=30.d0) !+ Specification and control of Amber's working precision !-------------BEGIN nmr.h ------------------------------------------------ ! ---Header file for the chemical shifts and NOESY intensity ! Because of the complexity of the storage requirements for ! these calculations, (and because one of the authors is very ! lazy,) storage allocation for these is not done in the ! LOCMEM routine, but rather arranged at compile time through ! the information given below. ! If you do not plan to make use of this section of the code, ! set the parameters below to small numbers to avoid allocating ! space unnecessarily. When/if you change your mind, reset the ! parameters and re-compile. ! Input parameters (things you have to set): ! MATOM = max # of atoms in the system ! MXR = max # residues ! MA = max # of protons in any sub-molecule ! MXTAU = max number of mixing times ! MXP = max number of input intensities (peaks) per mixing time ! MTOT = max. # total peaks, all mixing times, all sub-molecules ! MXVAR = max. # of "extra" dynamic variables ! MRING = max. # of rings for chemical shift calculation ! MSHF = max # of protons whose shifts are to be calculated ! MAXDIP = max # of residueal dipolar couplings ! MAXCSA = max # of residual csa measurements integer mring,mshf,mxvar,matom,ma,ma2,lst,mxr,mxtau,mxp, & mtot,maxdip,maxdipsets,maxcsa ! --- "standard" parameters for jobs like the test cases: parameter (mring=50) parameter (mshf=500) parameter (mxvar=10) parameter (matom=50000) parameter (ma=100) parameter (ma2=ma*ma) parameter (lst=(ma2+ma)/2) parameter (mxr=300) parameter (mxtau=5) parameter (mxp=100) parameter (mtot=500) parameter (maxdip=2000) parameter (maxcsa=200) parameter (maxdipsets=2) integer isubi,isubr parameter (isubi=8 + 3*ma + 2*matom + mxtau + 2*mxtau*mxp + mxp) parameter (isubr=3*ma + mxtau + 3*mxtau*mxp + 4) integer peakid(mxp) double precision tau(ma),pop(ma),popn(ma),emix(mxtau), & aexp(mxtau,mxp),awt(mxtau,mxp),arange(mxtau,mxp),oscale, & omega,taumet,taurot,invwt1,invwt2 integer nath,natmet,nummt,id2o,iroesy,ihet,nvect,ihsful,m2(ma), & inn(ma),ihyp(ma),ihyd(matom),inatom(matom),npeak(mxtau), & ihp(mxtau,mxp),jhp(mxtau,mxp) common /methylr/ tau,pop,popn,emix,aexp,awt,arange,oscale, & omega,taumet,taurot,invwt1,invwt2 common /methyli/ nath,natmet,nummt,id2o,iroesy,ihet,nvect,ihsful,m2, & inn,ihyp,ihyd,inatom,npeak,ihp,jhp,peakid ! Parameters for parallel broadcast integer BC_METHYLR,BC_METHYLI,BC_ALIGNR,BC_ALIGNI parameter(BC_METHYLR=3*ma+mxtau+3*mxtau*mxp+6) parameter(BC_METHYLI=8 + 3*ma + 2*matom + mxtau + 2*mxtau*mxp) parameter(BC_ALIGNR=5*maxdip + 1 + 6*maxdipsets) parameter(BC_ALIGNI=3*maxdip + 4) integer nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar,morse double precision scalm,pencut,ensave,tausw,ebdev,eadev,drjar common/nmr1/scalm,pencut,ensave,tausw,ebdev,eadev,drjar, & nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar,morse character(len=14) resat(matom) common/nmr2/ resat ! residual dipolar coupling (aka "alignment") restraint information: double precision dobsu(maxdip),dobsl(maxdip),dcut,gigj(maxdip), & dij(maxdip),dwt(maxdip), & s11(maxdipsets),s12(maxdipsets),s13(maxdipsets), & s22(maxdipsets),s23(maxdipsets),s33(maxdipsets) integer ndip,num_datasets,id(maxdip),jd(maxdip),dataset(maxdip),ifreeze, & ifreezes common/align/dobsu,dobsl,dcut,gigj,dij,dwt,s11,s12,s13,s22,s23,s33, & ndip,id,jd,dataset,num_datasets,ifreeze,ifreezes ! residual csa shift restraint information: double precision cobsu(maxcsa),cobsl(maxcsa),ccut,cwt(maxcsa), & sigma(3,3,maxcsa),field(maxcsa) integer ncsa,icsa(maxcsa),jcsa(maxcsa),kcsa(maxcsa),datasetc(maxcsa) common/csa/cobsu,cobsl,ccut,cwt,sigma,field,ncsa,icsa,jcsa,kcsa,datasetc integer mxvect parameter (mxvect=1) ! Common block containing variables relating to nmr restraints. integer intreq,irlreq,lnmr01,inmr02,iprr,iprw common/nmrstf/intreq,irlreq,lnmr01,inmr02,iprr,iprw integer ntot,ntota,ipmix(mtot),ntotb double precision calc(mtot),exper(mtot),calca(mtot),expera(mtot),calcb(mtot),experb(mtot) common/correl/ntot,ntota,calc,exper,calca,expera,calcb,experb,ipmix,ntotb double precision wnoesy,wshift,enoe,eshf,epcshf,ealign,ecsa common/wremar/wnoesy,wshift,enoe,eshf,epcshf,ealign,ecsa !-------------END nmr.h ------------------------------------------------ ! frames are centered on an atom atcenter that "owns" extra points ! numep are number of extra points attached to atcenter ! epframe are pointers to (at most 2) extra points attached to atcenter ! type is 1 if atom has at least 2 other bonds, pref to heavy atoms ! type is 2 if atom has only one other bond e.g. O6 in guanine ! first middle and third are atom nums defining frame ! loc_frame are local coords integer ifrtyp,iatcen,inumep,iepfr,ifrst,imid,ithrd integer leploc integer numfr,numextra integer numnb14,inb_14 common/fr_ptr/ifrtyp,iatcen,inumep,iepfr,ifrst,imid,ithrd, & leploc,inb_14,numfr,numextra,numnb14 !---------------------- les.h -------------------- !+ Specification and control of Amber's working precision ! parameters for 1: integer maxles,maxlestyp,maxlesadj parameter (maxles=500000) parameter (maxlestyp=100) parameter (maxlesadj=3000000) integer bc_lesr,bc_lesi parameter( bc_lesi=1+maxles*3+maxlesadj*2+1) parameter (bc_lesr=maxlestyp*maxlestyp+1) double precision lesfac(maxlestyp*maxlestyp),lfac ! for separate 1 and non-1 temperature coupling double precision ekinles0,temp0les,rndfles,sdfacles double precision scaltles,tempsules,ekeles,rsdles double precision ekmhles,ekphles common/lesr/lesfac,temp0les ! ileslst, jleslst and nlesadj are for PME and should maybe only be defined ! if the PME compile option is on integer ileslst(maxlesadj),jleslst(maxlesadj),nlesadj integer lestyp(maxles),nlesty,lestmp,cnum(maxles),subsp(maxles) ! this one is 1+maxles*3+maxlesadj*2+1 common/lesi/nlesty,lestyp,cnum,subsp,ileslst,jleslst,nlesadj ! some PME variables ! these are communicated in places other than parallel.f! but the sizes should ! not change double precision eeles,les_vir(3,3) common/ewlescomm/eeles,les_vir ! some partial REM variables double precision elesa,elesb,elesd,elesp common/eprem/elesa,elesb,elesd,elesp !---------------------- END les.h -------------------- integer nf integer i,nspsol,iok integer nhparm,idum,nttyp integer mbper,mgper,mdper,mbona,mtheta,mphia ! read but ignored character(len=80) fmt character(len=80) fmtin,ifmt,afmt,rfmt,type, line character(len=80), allocatable, dimension(:) :: ffdesc integer :: nlines, iscratch, ier ifmt = '(12I6)' afmt = '(20A4)' rfmt = '(5E16.8)' ! ----- READ THE MOLECULAR TOPOLOGY ----- nspsol = 0 ! ----- FORMATTED INPUT ----- fmtin = '(A80)' type = 'TITLE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmtin) title fmtin = ifmt type = 'POINTERS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) natom,ntypes,nbonh,mbona,ntheth,mtheta,nphih,mphia, & nhparm,nparm,nnb,nres,nbona,ntheta,nphia, & numbnd,numang,nptra,natyp,nphb,ifpert,nbper,ngper,ndper, & mbper,mgper,mdper,ifbox,nmxrs,ifcap,numextra,ncopy if (nparm /= 1) then write(6,*) ' *** THIS VERSION ONLY ACCEPTS TOPOLOGY FILES' write(6,*) ' THAT WERE CREATED BY ADDLES, WITH NPARM=1' write(6,*) ' USE A VERSION COMPILED WITHOUT -DLES ' call mexit(6,1) end if write(6,8118) natom,ntypes,nbonh,mbona,ntheth,mtheta,nphih,mphia, & nhparm,nparm,nnb,nres,nbona,ntheta,nphia,numbnd, & numang,nptra,natyp,nphb,ifbox,nmxrs,ifcap,numextra,ncopy 8118 format(t2, & 'NATOM = ',i7,' NTYPES = ',i7,' NBONH = ',i7,' MBONA = ',i7, & /' NTHETH = ',i7,' MTHETA = ',i7,' NPHIH = ',i7,' MPHIA = ',i7, & /' NHPARM = ',i7,' NPARM = ',i7,' NNB = ',i7,' NRES = ',i7, & /' NBONA = ',i7,' NTHETA = ',i7,' NPHIA = ',i7,' NUMBND = ',i7, & /' NUMANG = ',i7,' NPTRA = ',i7,' NATYP = ',i7,' NPHB = ',i7, & /' IFBOX = ',i7,' NMXRS = ',i7,' IFCAP = ',i7,' NEXTRA = ',i7 & ,/' NCOPY = ',i7/) ! --- make sure we do not exceed memory limits in commons --- nttyp = ntypes*(ntypes+1)/2 if (numbnd > MAX_BOND_TYPE .or. & numang > 900 .or. & nptra > 1200 .or. & nphb > 200 .or. & natyp > MAX_ATOM_TYPE .or. & nttyp > 1830 ) then write(6,'(/,5x,a)') 'rdparm: a parameter array overflowed' write(6,'(/,5x,a)') ' (e.g. the table of dihedral params)' write(6,8119) numbnd,MAX_BOND_TYPE, & numang,nptra, & natyp,MAX_ATOM_TYPE, & nphb,nttyp 8119 format( & ' NUMBND = ',i7,' max is',i5, & /' NUMANG = ',i7, ' max is 900', & /' NPTRA = ',i7, ' max is 1200', & /' NATYP = ',i7, ' max is',i5, & /' NPHB = ',i7, ' max is 200', & /' NTTYP = ',i7, ' max is 1830' & ) call mexit(6, 1) end if if(nbona /= mbona .or. ntheta /= mtheta .or. nphia /= mphia) then write(6,*) 'Sander no longer allows constraints in prmtop' write(6,*) '...must have nbona=mbona, ntheta=mtheta, nphi=mphi' call mexit(6,1) end if ! Write implicit solvent radius and screening info to mdout if (( igb /= 0 .and. (ifcap == 0 .or. ifcap == 5)).or.hybridgb>0) then fmtin = afmt type = 'RADIUS_SET' call nxtsec(nf, 6, 1,fmtin, type, fmt, iok) if (iok == 0) then ! Allow failure to support pre-AMBER9 prmtop read(nf,fmt) type ! Reuse type var to avoid declaring a new one write(6,'(A,A)') ' Implicit solvent radii are ',type end if if ( igb == 7 ) then write(6,'(A)') ' Replacing prmtop screening parameters with GBn (igb=7) values' end if end if ! --- Read the force field information from the prmtop if available fmtin = '(i,a)' type = 'FORCE_FIELD_TYPE' call nxtsec(nf, 6, 1,fmtin, type, fmt, iok) if(iok == 0) then ! We found a force field description. Should be 1 or more lines ! of text which will be echoed back to the user. In some case ! we also take some decisions based on this. write (6,'(a)') '| Force field information read from topology file: ' ! format should be (nlines, text) where we read the first lines nlines ! to know how many more lines to read - we ignore nlines on subsequent ! lines but they still need to be there. ! e.g. ! 2 PARM99 ! 2 GAFF read(nf,fmt) nlines,line allocate (ffdesc(nlines), stat = ier) if(.not.(ier==0))call Aass('ier==0',"rdparm.f",145) ffdesc(1) = line write(6,'(a,a)') '| ',ffdesc(1) do i = 2, nlines read(nf,fmt) iscratch,ffdesc(i) write(6,'(a,a)') '| ',ffdesc(i) end do !Test to see if any special force fields are in use. !1) Check for CHARMM charmm = .false. do i = 1, nlines if (index(ffdesc(i),'CHARMM') /= 0) then write(6,'("|")') write(6,'(a,a)') '| CHARMM force field in use. ' charmm = .true. !AMBER 10 does not support CHARMM prmtop files quit !with an error message. call sander_bomb('rdparm','CHARMM prmtop file detected.', & 'AMBER 10 does not support CHARMM prmtop files.') !test to see if scee and scnb are set to 1.0 if (scee /= 1.0d0) then write(6,'(a)') ' WARNING: CHARMM force field is in use but scee /= 1.0.' write(6,'(a,f4.2)') ' scee = ',scee write(6,'(a)') ' results will be INCORRECT!' if(.not.(scee==1.0d0))call Aass('scee==1.0d0',"rdparm.f",173) end if if (scnb /= 1.0d0) then write(6,'(a)') ' WARNING: CHARMM force field is in use but scnb /= 1.0.' write(6,'(a,f4.2)') ' scnb = ',scnb write(6,'(a)') ' results will be INCORRECT!' if(.not.(scnb==1.0d0))call Aass('scnb==1.0d0',"rdparm.f",179) end if end if end do if(charmm)then type = 'CHARMM_NUM_IMPROPERS' call nxtsec(nf, 6, 1,fmtin, type, fmt, iok) read(nf,fmt) nimphi write(0,*)"Got number of impropers: ",nimphi allocate(im(nimphi),jm(nimphi),km(nimphi),lm(nimphi),imp(nimphi), & stat=ier) if(.not.(ier==0))call Aass('ier==0',"rdparm.f",191) endif !End Test deallocate (ffdesc, stat = ier) if(.not.(ier==0))call Aass('ier==0',"rdparm.f",197) end if return end subroutine rdparm1 !======================================================================= !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine rdparm2 here] subroutine rdparm2(x,ix,ih,ipairs,nf) use parms use nblist, only: a,b,c use amoeba_mdin,only: iamoeba use pimd_vars, only: ipimd, dmdlm, itimass implicit none integer nf double precision x(*) integer ix(*),ipairs(*) character(len=4) ih(*) !+ Specification and control of Amber's working precision !-------------BEGIN md.h ------------------------------------------------ integer BC_MDI ! size in integers of common block mdi integer BC_MDR ! size in Reals 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,alpb,iyammp, &!41 gbsa,vrand,iwrap,nrespa,irespa,nrespai,icfe, &!48 rbornstat,ivcap,iconstreff, &!51 neb,vv,tmode,ipol,iesp,ievb,nodeid,num_noshake, &!59 idecomp,icnstph,ntcnstph,maxdup,numexchg,repcrd,numwatkeep,hybridgb, &!67 ibgion,ienion,profile_mpi !70 parameter (BC_MDI=70) 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,alpb,iyammp,gbsa,vrand,numexchg,repcrd,numwatkeep,hybridgb, & iwrap,nrespa,irespa,nrespai,icfe,rbornstat, & ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup,neb,vv, & tmode,ipol,iesp,ievb,nodeid,num_noshake,ibgion,ienion, profile_mpi ! ... floats: double precision t,dt,temp0,tautp,pres0,comp,taup,temp,tempi, & !9 tol,taur,dx0,drms,vlimit,rbtarg(9),tmass,tmassinv, & !25 kappa,offset,surften,gamma_ln,extdiel,intdiel,rdt, & !32 gbalpha,gbbeta,gbgamma,cut_inner,clambda,saltcon, & !38 solvph,rgbmax,fsmax,restraint_wt, & !42 skmin,skmax,vfac,gbneckscale,v11,v12,v22,kevb,evbt,Arad !52 parameter (BC_MDR=52) 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, & gbalpha,gbbeta,gbgamma,cut_inner,clambda,saltcon, & solvph,rgbmax,fsmax,restraint_wt,skmin,skmax,vfac,gbneckscale, & v11,v12,v22,kevb,evbt,Arad ! ... strings: character(len=4) iwtnm,iowtnm,ihwtnm character(len=256) restraintmask,bellymask,tgtfitmask,tgtrmsmask,noshakemask,crgmask common/mds/ restraintmask,bellymask,tgtfitmask,tgtrmsmask,noshakemask,crgmask, & iwtnm,iowtnm,ihwtnm(2) !-------------END md.h ------------------------------------------------ ! --- following should not need to be modified unless you are adding ! more variables to the "locmem" style of memory management ! 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,nparm, & natc,nattgtfit,nattgtrms,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,itgtfitgp,itgtrmsgp,i64,i65,i68, & i70,i72,i74,i76,i78,i79,i80,i82,i84,i86, & icpstinf,icpresst,icptrsct, icpptcnt, & l15,lwinv,lpol,lcrd,lforce,l36,lvel,lvel2,l45,l50, & lcrdr,l60,l65,lmass,l75,l80,l85,l90,l95,l96,l97,l98,l99,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, & iifstwt,iifstwr,nrealb,nintb,nholb,npairb,lastr,lasti,lasth, & lastpr,nbper,ngper,ndper,ifpert,ncopy, & imask1,imask2,numadjst,mxadjmsk,icphidx,icptpair,lfsg,lvsg,noshake ! 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,nparm , & !14 natc ,nattgtfit,nattgtrms,ibelly ,natbel,ishake,nmxrs , & !21 mxsub ,natyp ,npdec ,i02 ,i04 ,i06 ,i08 ,i10 , & !29 iibh ,ijbh ,iicbh ,iiba ,ijba ,iicba , & !35 i24 ,i26 ,i28 ,i30 ,i32 ,i34 ,i36 ,i38 ,i40 , & !44 i42 ,i44 ,i46 ,i48 ,i50 ,i52 ,i54 ,i56 ,i58 ,ibellygp,& !54 icnstrgp,itgtfitgp,itgtrmsgp,i64 ,i65 ,i68 , & !60 i70 ,i72 ,i74 ,i76 ,i78 ,i79 ,i80 ,i82 , & !68 i84 ,i86 , & !70 icpstinf,icpresst ,icptrsct ,icpptcnt, & !74 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 ,l99 ,lfrctmp , & !98 l105 ,l110 ,l115 ,l120 ,l125 ,l130 ,l135 ,l140 ,l145 ,l150 , & !108 l165 ,l170 ,l175 ,l180 ,l185 ,l186 ,l187 ,l188 ,l189 ,l190 , & !118 lcpcrg ,lcpene , & !120 m02 ,m04 ,m06 ,m08 ,m10 ,m12 ,m14 ,m16 ,m18 ,i01 , & !130 iifstwt ,iifstwr ,nrealb ,nintb ,nholb ,npairb,lastr ,lasti ,lasth , & !139 lastpr ,nbper ,ngper ,ndper ,ifpert, ncopy , & !145 imask1 ,imask2 ,numadjst ,mxadjmsk,icphidx,icptpair,lfsg,lvsg,noshake !154 ! 154 is the size of the MEMORY common block !+ Specification and control of Amber's Input/Output ! File names character(len=4096) groupbuffer character(len=256) mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmrf, mincor, & vecs, radii, freqe,redir(9),rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt, evbin, evbout, mmtsb_setup_file,pimdout, & inptraj character owrite, facc common /files/ groupbuffer, mdin, mdout, inpcrd, parm, restrt, & refc, mdvel, mden, mdcrd, mdinfo, nmrf, mincor, & vecs, radii, freqe, owrite, facc,rstdip,mddip,inpdip,groups,gpes, & cpin, cpout, cprestrt, evbin, evbout, mmtsb_setup_file,pimdout, & inptraj ! put this in a separate common block to stop the compiler from ! complaining about misalignment integer numgroup, nslice common/nmgrp/ numgroup, nslice ! File units ! An I/O Unit resource manager does not exist. integer MDCRD_UNIT integer INPTRAJ_UNIT integer MDEN_UNIT integer MDINFO_UNIT integer MDVEL_UNIT parameter ( MDINFO_UNIT = 7 ) parameter ( MDCRD_UNIT = 12 ) parameter ( INPTRAJ_UNIT = 24 ) 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 ) !! !! EVB I/O unit !! integer, parameter :: evb_unit = 75 integer, parameter :: schlegel_unit = 80 !! FULL PIMD I/O unit integer, parameter :: pimd_unit = 277 ! File related controls and options character(len=80) title,title1 common/runhed/ title, title1 logical mdin_ewald,mdin_pb,mdin_amoeba common/mdin_flags/mdin_ewald,mdin_pb,mdin_amoeba integer BC_HULP ! size in integers of common HULP parameter ( BC_HULP = 9 ) integer ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave common/hulp/ntpr,ntwr,ntwx,ntwv,ntwe,ntpp,ioutfm,ntwprt,ntave ! 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(9) common/nmrrdr/redir,iredir ! common block sizes: ! ... floats: !+ Specification and control of Amber's working precision double precision box,cut,scnb,scee,dielc,rad,wel,radhb,welhb, & cutcap,xcap,ycap,zcap,fcap, & xlorth,ylorth,zlorth,xorth,yorth,zorth,forth, & rwell,xbox0,ybox0,zbox0 common/boxr/box(3),cut,scnb,scee,dielc,xbox0,ybox0,zbox0, & cutcap,xcap,ycap,zcap,fcap, & xlorth,ylorth,zlorth,xorth,yorth,zorth,forth, & rwell,rad(100),wel(100),radhb(200),welhb(200) ! ... integers: integer ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp common/boxi/ntb,ifbox,numpk,nbit,ifcap,natcap,isftrp double precision extraboxdim parameter (extraboxdim=30.d0) !+ Specification and control of Amber's working precision !-------------BEGIN nmr.h ------------------------------------------------ ! ---Header file for the chemical shifts and NOESY intensity ! Because of the complexity of the storage requirements for ! these calculations, (and because one of the authors is very ! lazy,) storage allocation for these is not done in the ! LOCMEM routine, but rather arranged at compile time through ! the information given below. ! If you do not plan to make use of this section of the code, ! set the parameters below to small numbers to avoid allocating ! space unnecessarily. When/if you change your mind, reset the ! parameters and re-compile. ! Input parameters (things you have to set): ! MATOM = max # of atoms in the system ! MXR = max # residues ! MA = max # of protons in any sub-molecule ! MXTAU = max number of mixing times ! MXP = max number of input intensities (peaks) per mixing time ! MTOT = max. # total peaks, all mixing times, all sub-molecules ! MXVAR = max. # of "extra" dynamic variables ! MRING = max. # of rings for chemical shift calculation ! MSHF = max # of protons whose shifts are to be calculated ! MAXDIP = max # of residueal dipolar couplings ! MAXCSA = max # of residual csa measurements integer mring,mshf,mxvar,matom,ma,ma2,lst,mxr,mxtau,mxp, & mtot,maxdip,maxdipsets,maxcsa ! --- "standard" parameters for jobs like the test cases: parameter (mring=50) parameter (mshf=500) parameter (mxvar=10) parameter (matom=50000) parameter (ma=100) parameter (ma2=ma*ma) parameter (lst=(ma2+ma)/2) parameter (mxr=300) parameter (mxtau=5) parameter (mxp=100) parameter (mtot=500) parameter (maxdip=2000) parameter (maxcsa=200) parameter (maxdipsets=2) integer isubi,isubr parameter (isubi=8 + 3*ma + 2*matom + mxtau + 2*mxtau*mxp + mxp) parameter (isubr=3*ma + mxtau + 3*mxtau*mxp + 4) integer peakid(mxp) double precision tau(ma),pop(ma),popn(ma),emix(mxtau), & aexp(mxtau,mxp),awt(mxtau,mxp),arange(mxtau,mxp),oscale, & omega,taumet,taurot,invwt1,invwt2 integer nath,natmet,nummt,id2o,iroesy,ihet,nvect,ihsful,m2(ma), & inn(ma),ihyp(ma),ihyd(matom),inatom(matom),npeak(mxtau), & ihp(mxtau,mxp),jhp(mxtau,mxp) common /methylr/ tau,pop,popn,emix,aexp,awt,arange,oscale, & omega,taumet,taurot,invwt1,invwt2 common /methyli/ nath,natmet,nummt,id2o,iroesy,ihet,nvect,ihsful,m2, & inn,ihyp,ihyd,inatom,npeak,ihp,jhp,peakid ! Parameters for parallel broadcast integer BC_METHYLR,BC_METHYLI,BC_ALIGNR,BC_ALIGNI parameter(BC_METHYLR=3*ma+mxtau+3*mxtau*mxp+6) parameter(BC_METHYLI=8 + 3*ma + 2*matom + mxtau + 2*mxtau*mxp) parameter(BC_ALIGNR=5*maxdip + 1 + 6*maxdipsets) parameter(BC_ALIGNI=3*maxdip + 4) integer nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar,morse double precision scalm,pencut,ensave,tausw,ebdev,eadev,drjar common/nmr1/scalm,pencut,ensave,tausw,ebdev,eadev,drjar, & nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar,morse character(len=14) resat(matom) common/nmr2/ resat ! residual dipolar coupling (aka "alignment") restraint information: double precision dobsu(maxdip),dobsl(maxdip),dcut,gigj(maxdip), & dij(maxdip),dwt(maxdip), & s11(maxdipsets),s12(maxdipsets),s13(maxdipsets), & s22(maxdipsets),s23(maxdipsets),s33(maxdipsets) integer ndip,num_datasets,id(maxdip),jd(maxdip),dataset(maxdip),ifreeze, & ifreezes common/align/dobsu,dobsl,dcut,gigj,dij,dwt,s11,s12,s13,s22,s23,s33, & ndip,id,jd,dataset,num_datasets,ifreeze,ifreezes ! residual csa shift restraint information: double precision cobsu(maxcsa),cobsl(maxcsa),ccut,cwt(maxcsa), & sigma(3,3,maxcsa),field(maxcsa) integer ncsa,icsa(maxcsa),jcsa(maxcsa),kcsa(maxcsa),datasetc(maxcsa) common/csa/cobsu,cobsl,ccut,cwt,sigma,field,ncsa,icsa,jcsa,kcsa,datasetc integer mxvect parameter (mxvect=1) ! Common block containing variables relating to nmr restraints. integer intreq,irlreq,lnmr01,inmr02,iprr,iprw common/nmrstf/intreq,irlreq,lnmr01,inmr02,iprr,iprw integer ntot,ntota,ipmix(mtot),ntotb double precision calc(mtot),exper(mtot),calca(mtot),expera(mtot),calcb(mtot),experb(mtot) common/correl/ntot,ntota,calc,exper,calca,expera,calcb,experb,ipmix,ntotb double precision wnoesy,wshift,enoe,eshf,epcshf,ealign,ecsa common/wremar/wnoesy,wshift,enoe,eshf,epcshf,ealign,ecsa !-------------END nmr.h ------------------------------------------------ !---------------------- les.h -------------------- !+ Specification and control of Amber's working precision ! parameters for 1: integer maxles,maxlestyp,maxlesadj parameter (maxles=500000) parameter (maxlestyp=100) parameter (maxlesadj=3000000) integer bc_lesr,bc_lesi parameter( bc_lesi=1+maxles*3+maxlesadj*2+1) parameter (bc_lesr=maxlestyp*maxlestyp+1) double precision lesfac(maxlestyp*maxlestyp),lfac ! for separate 1 and non-1 temperature coupling double precision ekinles0,temp0les,rndfles,sdfacles double precision scaltles,tempsules,ekeles,rsdles double precision ekmhles,ekphles common/lesr/lesfac,temp0les ! ileslst, jleslst and nlesadj are for PME and should maybe only be defined ! if the PME compile option is on integer ileslst(maxlesadj),jleslst(maxlesadj),nlesadj integer lestyp(maxles),nlesty,lestmp,cnum(maxles),subsp(maxles) ! this one is 1+maxles*3+maxlesadj*2+1 common/lesi/nlesty,lestyp,cnum,subsp,ileslst,jleslst,nlesadj ! some PME variables ! these are communicated in places other than parallel.f! but the sizes should ! not change double precision eeles,les_vir(3,3) common/ewlescomm/eeles,les_vir ! some partial REM variables double precision elesa,elesb,elesd,elesp common/eprem/elesa,elesb,elesd,elesp !---------------------- END les.h -------------------- integer iexcl,numex,k1,j1 integer nttyp,ntype,i,iok integer i1,i2,i3,i4,j,jj,k,l,n,nn integer iptres,nspsol,natsm,idum,ip14 integer l_ib,l_jb,l_bt integer l_it,l_jt,l_kt,l_tt integer l_id,l_jd,l_kd,l_ld,l_dt integer bp,ibp,jbp,btp double precision dumd,oldbeta,duma,dumb,dumc character(len=4) dumchar integer dumint double precision dumfloat double precision massdiff ! = mass[perturbed] - mass[original] for TI w.r.t. mass. integer irotat( natom ) ! dummy: disappears when leaving subroutine integer allocate_err integer ierr ! Allocation status. character(len=80) fmt character(len=80) fmtin,ifmt,afmt,rfmt,type ifmt = '(12I6)' afmt = '(20A4)' rfmt = '(5E16.8)' nttyp = ntypes*(ntypes+1)/2 ntype = ntypes*ntypes ! ----- READ THE SYMBOLS AND THE CHARGES AND THE MASSES ----- fmtin = afmt type = 'ATOM_NAME' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ih(m04+i-1),i = 1,natom) fmtin = rfmt type = 'CHARGE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (x(l15+i-1),i = 1,natom) fmtin = rfmt type = 'MASS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (x(lwinv+i-1),i = 1,natom) fmtin = ifmt type = 'ATOM_TYPE_INDEX' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i04+i-1),i = 1,natom) fmtin = ifmt type = 'NUMBER_EXCLUDED_ATOMS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i08-1),i = 1,natom) fmtin = ifmt type = 'NONBONDED_PARM_INDEX' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i06-1),i = 1,ntype) fmtin = afmt type = 'RESIDUE_LABEL' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ih(i+m02-1),i=1,nres) fmtin = ifmt type = 'RESIDUE_POINTER' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i02-1),i=1,nres) ix(i02+nres) = natom+1 ! ----- READ THE PARAMETERS ----- fmtin = rfmt type = 'BOND_FORCE_CONSTANT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (rk(i), i = 1,numbnd) fmtin = rfmt type = 'BOND_EQUIL_VALUE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (req(i), i = 1,numbnd) fmtin = rfmt type = 'ANGLE_FORCE_CONSTANT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (tk(i), i = 1,numang) fmtin = rfmt type = 'ANGLE_EQUIL_VALUE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (teq(i), i = 1,numang) fmtin = rfmt type = 'DIHEDRAL_FORCE_CONSTANT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (pk(i), i = 1,nptra) fmtin = rfmt type = 'DIHEDRAL_PERIODICITY' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (pn(i), i = 1,nptra) fmtin = rfmt type = 'DIHEDRAL_PHASE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (phase(i), i = 1,nptra) fmtin = rfmt type = 'SOLTY' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (solty(i), i = 1,natyp) fmtin = rfmt type = 'LENNARD_JONES_ACOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (cn1(i), i = 1,nttyp) fmtin = rfmt type = 'LENNARD_JONES_BCOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (cn2(i), i = 1,nttyp) ! ----- READ THE BONDING INFORMATIONS ----- fmtin = ifmt type = 'BONDS_INC_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+iibh-1),ix(i+ijbh-1),ix(i+iicbh-1), & i = 1,nbonh) fmtin = ifmt type = 'BONDS_WITHOUT_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt)(ix(i+iiba-1),ix(i+ijba-1),ix(i+iicba-1),i = 1,nbona) fmtin = ifmt type = 'ANGLES_INC_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i24-1),ix(i+i26-1),ix(i+i28-1),ix(i+i30-1), & i = 1,ntheth) fmtin = ifmt type = 'ANGLES_WITHOUT_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i32-1),ix(i+i34-1),ix(i+i36-1),ix(i+i38-1), & i = 1,ntheta) fmtin = ifmt type = 'DIHEDRALS_INC_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i40-1),ix(i+i42-1),ix(i+i44-1),ix(i+i46-1), & ix(i+i48-1),i = 1,nphih) fmtin = ifmt type = 'DIHEDRALS_WITHOUT_HYDROGEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i50-1),ix(i+i52-1),ix(i+i54-1),ix(i+i56-1), & ix(i+i58-1),i = 1,nphia) fmtin = ifmt type = 'EXCLUDED_ATOMS_LIST' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i10-1),i=1,nnb) ! ----- READ THE H-BOND PARAMETERS ----- fmtin = rfmt type = 'HBOND_ACOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (asol(i),i=1,nphb) do i=1,nphb if( asol(i) /= 0.d0 ) then write(6,*) 'Found a non-zero 10-12 coefficient, but source', & ' was not compiled with -DHAS_10_12.' write(6,*) 'If you are using a pre-1994 force field, you', & ' will need to re-compile with this flag.' call mexit(6,1) end if end do fmtin = rfmt type = 'HBOND_BCOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (bsol(i),i=1,nphb) fmtin = rfmt type = 'HBCUT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (hbcut(i),i=1,nphb) ! ----- READ ISYMBL,ITREE,JOIN AND IROTAT ARRAYS ----- fmtin = afmt type = 'AMBER_ATOM_TYPE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ih(i+m06-1),i=1,natom) fmtin = afmt type = 'TREE_CHAIN_CLASSIFICATION' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ih(i+m08-1),i=1,natom) fmtin = ifmt type = 'JOIN_ARRAY' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i64-1),i=1,natom) fmtin = ifmt type = 'IROTAT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (irotat(i), i=1,natom) ! ----- READ THE BOUNDARY CONDITION STUFF ----- nspm = 1 ix(i70) = natom if (ifbox > 0) then fmtin = ifmt type = 'SOLVENT_POINTERS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) iptres,nspm,nspsol fmtin = ifmt type = 'ATOMS_PER_MOLECULE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (ix(i+i70-1),i=1,nspm) fmtin = rfmt type = 'BOX_DIMENSIONS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) oldbeta,duma,dumb,dumc ! ---(above values are still read, for backward compatibility, but ! ignored. Box info must come from the coord. file or from ! the &ewald namelist of the input file) if( igb /= 0 .or. ntb == 0 )then box(1)=0.0d0 box(2)=0.0d0 box(3)=0.0d0 else box(1)=a box(2)=b box(3)=c end if end if ! (ifbox > 0) ! ----- LOAD THE CAP INFORMATION IF NEEDED ----- if(ifcap == 1) then fmtin = '(I6)' type = 'CAP_INFO' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) natcap fmtin = '(4E16.8)' type = 'CAP_INFO2' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) cutcap,xcap,ycap,zcap end if if( igb /= 0 .and. ifcap == 0 .and. iok == -1 ) then write(0,*) 'GB/PB calculations now require a new-style prmtop file' write(6,*) 'GB/PB calculations now require a new-style prmtop file' call mexit(6,1) end if if (( igb /= 0 .and. (ifcap == 0 .or. ifcap == 5)).or.hybridgb>0) then fmtin = rfmt type = 'RADII' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (x(l97+i-1),i=1,natom) type = 'SCREEN' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (x(l96+i-1),i=1,natom) end if if (ipol > 0) then fmtin = rfmt type = 'POLARIZABILITY' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (x(lpol+i-1),i=1,natom) end if ! ----- READ THE PERTURBED MASSES IF NEEDED ----- if (itimass > 0) then ! JVAN: dmdlm must be allocated here, not in pimd_init. allocate( dmdlm(1:natom),stat=ierr ) if(.not.( ierr == 0 ))call Aass(' ierr == 0 ',"rdparm.f",523) fmtin = rfmt type = 'TI_MASS' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (dmdlm(i) ,i = 1,natom) do i=1,natom massdiff = (dmdlm(i) - x(lwinv+i-1)) x(lwinv+i-1) = x(lwinv+i-1) + clambda * massdiff dmdlm(i) = massdiff/x(lwinv+i-1) end do end if if (nparm == 1.and.iamoeba.eq.0) then fmtin = ifmt type = 'LES_NTYP' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read (nf,fmt) nlesty lestmp=nlesty*nlesty ! check the array sizes to make sure we do not overflow ! 1 types if (nlesty > maxlestyp) then write (6,*) 'Exceeded MAXLESTYP',nlesty stop end if ! 1 atoms if (natom > maxles) then write (6,*) 'Exceeded MAXLES',natom stop end if fmtin = ifmt type = 'LES_TYPE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read (nf,fmt) (lestyp(i),i=1,natom) fmtin = rfmt type = 'LES_FAC' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read (nf,fmt) (lesfac(i),i=1,lestmp) fmtin = ifmt type = 'LES_CNUM' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read (nf,fmt) (cnum(i), i=1,natom) fmtin = ifmt type = 'LES_ID' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read (nf,fmt) (subsp(i), i=1,natom) write (6,*) 'LES parameters were found' ! now create the list of atoms that have non-unitary scaling factors. ! this will be used in the Ewald calculation to correct for the ! lack of use of the intra-copy scaling factor in the charge grid. ! all of these pairs will need correction. The list will not change and ! is therefore only calculated once (here). nlesadj=0 iexcl=0 ! pairs are listed in two arrays, for i and j, rather than using ! a set of pointers like the nonbond and exclusion lists. This is ! since many atoms will not have any correction partners (since ! they are not in 1). if( ipimd.eq.0 ) then ! pimd are not going to use nb_adjust_les, so we do not need to generate ! les adjust list ! do 6500 k=1,natom lestmp=nlesty*(lestyp(k)-1) ! write (6,*) 'atom1 : ',k,lestmp ! need to sum all f the number of exclusions even if non-1 atoms. ! see below. numex=ix(k+i08-1) do 6510 j=k+1,natom lfac=lesfac(lestmp+lestyp(j)) ! check for non-zero scaling factor (meaning a correction will ! be required) if (abs(lfac-1.0d0) > 0.01) then ! check to make sure these aren't excluded atoms (since then ! no correction is wanted) ! FORMAT(12I6) (NATEX(i), i=1,NEXT) ! the excluded atom list. To get the excluded list for atom ! "i" you need to traverse the NUMEX list, adding up all ! the previous NUMEX values, since NUMEX(i) holds the number ! of excluded atoms for atom "i", not the index into the ! NATEX list. Let IEXCL = SUM(NUMEX(j), j=1,i-1), then ! excluded atoms are NATEX(IEXCL)+1 to NATEX(IEXCL+NUMEX(i)). do k1=iexcl+1,iexcl+numex ! get exclusion j1=ix(k1+i10-1) ! check against atom j if (j1 == j) then ! excluded, get next atom j ! write (6,*) 'Exclusion list match' goto 6510 end if ! check next entry end do ! if we arrived here, the atom was not in the exclusion list ! so this pair will need correction ! (should add boundary checking for variables here) if (nlesadj == maxlesadj) then write (6,*) 'EXCEEDED MAXLESADJ!' stop end if nlesadj=nlesadj+1 ileslst(nlesadj)=k jleslst(nlesadj)=j end if ! next j 6510 continue ! increment the exclusion list pointer for atom i iexcl=iexcl+numex ! next i 6500 continue end if !(ipimd == 0 ) write (6,6520) nlesadj 6520 format (1x,i7,' LES atom pairs require adjustment') ! end creation of 1 adjustment list and reading 1 info end if ! (nparm==1) ! ----- CALCULATE INVERSE, TOTAL MASSES ----- ! -- save the masses for removal of mass weighted velocity, ! leaving the inverse masses in the legacy, Lwinv area tmass = 0.0d0 ! -- index over molecules j = l75-1 jj = i70-1 ! -- index over mass->invmass k = lwinv-1 ! -- index over saved mass l = lmass-1 do n = 1,nspm j = j + 1 jj = jj + 1 x(j) = 0.0d0 natsm = ix(jj) do nn = 1,natsm k = k+1 l = l+1 ! -- sum molecule x(j) = x(j) + x(k) ! -- save mass in "new" Lmass area x(l) = x(k) ! -- make inverse in "old" Lwinv area if( x(k) /= 0.d0 ) x(k) = 1.0d0 / x(k) end do tmass = tmass + x(j) end do tmassinv = 1.0d0 / tmass ! ----- SCALE THE CHARGES IF DIELC.GT.1.0E0 ----- if (dielc /= 1.0e0 .and. igb == 0) then dumd = sqrt(dielc) do i = 1,natom x(i+l15-1) = x(i+l15-1)/dumd end do end if ! ----- INVERT THE HBCUT ARRAY ----- do i = 1,nphb if(hbcut(i) <= 0.001e0) hbcut(i) = 1.0d-10 hbcut(i) = 1.0e0/hbcut(i) end do ! ----- duplicate dihedral pointers for vector ephi ----- call dihdup(nphih,ix(i40),ix(i42),ix(i44),ix(i46),ix(i48),pn) call dihdup(nphia,ix(i50),ix(i52),ix(i54),ix(i56),ix(i58),pn) ! --- pre-calculate some parameters for vector ephi --- call dihpar(nptra,pk,pn,phase,gamc,gams,ipn,fmn) if (charmm) then ! ---read in 1-4 parameters ! allocate(cn114(nttyp),cn214(nttyp),stat=allocate_err) ! if(allocate_err /= 0) & ! call sander_bomb('rdparm2', & ! 'cannot allocate urey-bradley arrays',' BOMBING') fmtin = rfmt type = 'LENNARD_JONES_14_ACOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (cn114(i), i = 1,nttyp) fmtin = rfmt type = 'LENNARD_JONES_14_BCOEF' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (cn214(i), i = 1,nttyp) ! ! --- read in urey-bradley parameters: ! allocate(rub(numang),rkub(numang),stat=allocate_err) ! if(allocate_err /= 0) & ! call sander_bomb('rdparm2', & ! 'cannot allocate urey-bradley arrays',' BOMBING') fmtin = rfmt type = 'UREY_BRADLEY_EQUIL_VALUE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (rub(i), i = 1,numang) fmtin = rfmt type = 'UREY_BRADLEY_FORCE_CONSTANT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (rkub(i), i = 1,numang) fmtin = '(i8)' type = 'CHARMM_NUM_IMPR_TYPES' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) nimprtyp write(0,*)"Got number of improper types: ",nimprtyp fmtin = rfmt type = 'CHARMM_IMPROPER_FORCE_CONSTANT' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (pk_impr(i), i = 1,nimprtyp) write(0,fmt)(pk_impr(i), i = 1,nimprtyp) type = 'CHARMM_IMPROPER_PHASE' call nxtsec(nf, 6, 0,fmtin, type, fmt, iok) read(nf,fmt) (phase_impr(i), i = 1,nimprtyp) write(0,fmt)(phase_impr(i), i = 1,nimprtyp) end if return end subroutine rdparm2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine istuff here] subroutine istuff(i,j,iarray,k) ! routine to correctly load a strange shaped 2 dim matrix dimension iarray(15,*) iarray(i,j) = k return end subroutine istuff