C C $Id: mcsetup.F,v 1.13 1998/07/16 16:40:05 jjv5 Exp $ C C------------------------------------------------------------------------ subroutine mcsetup(wchar,bchar,newsub,wrtstuff,ok,rok, & ecrt,dcrt,ires,movedres,imc,igtyp,nstp,iter,ierror) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" character wchar, bchar logical newsub, wrtstuff, ok, rok dimension ires(maxres) C strip dr, dangle, seed, betapme, dbox and mc out of keywords C (since these change during mc-run and need to be filled out by C wrtdiv) call stripkeymc wchar = 'W' bchar = ' ' iter = 1 nscf = 0 ntry = 1 naccpt = 1 imc = 1 newsub = .true. wrtstuff = .true. wrtedmx = .true. ok = .true. rok = .true. pterm = 0.0 enth = 0.0 C set the convergence criteria ecrt = eecrt dcrt = dencrt C open some files if(myid.eq.0) then call opnfil(7,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(7),'(A80,"; NODE ",I4)') title, myid call opnfil(8,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(8),'(A80,"; NODE ",I4)') title, myid if (gradient) then call opnfil(10,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(10),'(A80,"; NODE ",I4)') title, myid if (xyzspc) then igtyp = 2 C-RDC write(iunit(10),'(" CARTESIAN COORDINATE GRADIENT:")') else igtyp = 1 C-RDC write(iunit(10),'(" INTERNAL COORDINATE GRADIENT:")') endif C-RDC write(iunit(10),'(" PARAMETER ATOM ELEMENTAL PARAMETER", C-RDC & 6X,"PARAMETER"/4X,"NO.",6X,"NO.",3X,"SYMBOL",6X, C-RDC & "TYPE",11X,"VALUE",13X,"GRADIENT"/)') endif else C call opnpfil(11,ierror) C if (ierror.ne.0) goto 1000 call opnpfil(17,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(17),'(A80,"; NODE ",I4)') title, myid call opnpfil(18,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(18),'(A80,"; NODE ",I4)') title, myid if (gradient) then call opnpfil(20,ierror) if (ierror.ne.0) goto 1000 C-RDC write(iunit(20),'(A80,"; NODE ",I4)') title, myid if (xyzspc) then igtyp = 2 C-RDC write(iunit(20),'(" CARTESIAN COORDINATE GRADIENT:")') else igtyp = 1 C-RDC write(iunit(20),'(" INTERNAL COORDINATE GRADIENT:")') endif C-RDC write(iunit(20),'(" PARAMETER ATOM ELEMENTAL PARAMETER", C-RDC & 6X,"PARAMETER"/4X,"NO.",6X,"NO.",3X,"SYMBOL",6X, C-RDC & "TYPE",11X,"VALUE",13X,"GRADIENT"/)') endif endif C initialize coordinates etc. call setup(ierror) if (ierror.ne.0) goto 1000 C get geometry centers of the residues and the boxdimensions call setbox(ierror) if (ierror.ne.0) goto 1000 C setup for gradient/virial call setopt C setup for PME if (pme) then call pme_setup call pme_calcb call pme_calctheta C-RDC write(iout,'(/" A CLASSICAL SMOOTH PARTICLE MESH EWALD ", C-RDC & "SUM (PME) IS USED TO", C-RDC & /" CALCULATE THE LONG RANGE ELECTROSTATIC INTERACTION.", C-RDC & //" PME-PARAMETERS:",/" --------------")') C-RDC write(iout,'(" NUMBER OF GRIDS IN THE X-DIRECTION: ", C-RDC & I4,/" NUMBER OF GRIDS IN THE Y-DIRECTION: ", I4, C-RDC & /" NUMBER OF GRIDS IN THE Z-DIRECTION: ", I4, C-RDC & /" ORDER OF CARDINAL B-SPLINE INTERPOLATION: ", I4, C-RDC & /" EXPONENTIAL BETA-FACTOR: ", F10.6)') C-RDC & k1pme,k2pme,k3pme,nspline,betapme if (cm1) then C-RDC write(iout,'(/" PME WILL USE CM1 CHARGES")') elseif (cm2) then C-RDC write(iout,'(/" PME WILL USE CM2 CHARGES")') else C-RDC write(iout,'(/" PME WILL USE MULLIKEN CHARGES")') endif if (recipintr) then C-RDC write(iout,'(/" ALL INTRAMOLECULAR CONTRIBUTIONS", C-RDC & " TO THE GRADIENT ARE NEGLECTED,", C-RDC & /" EXCEPT FOR THE INTRAMOLECULAR CONTRIBUTIONS", C-RDC & " OF THE RECIPROCAL ENERGY")') else C-RDC write(iout,'(/" ALL INTRAMOLECULAR CONTRIBUTIONS", C-RDC & " TO THE GRADIENT ARE NEGLECTED")') endif else C-RDC write(iout,'(/" ALL INTRAMOLECULAR CONTRIBUTIONS", C-RDC & " TO THE GRADIENT ARE NEGLECTED")') endif if (wrtmc) then C-RDC write(iscr,'(/"****************************", C-RDC & "***********************", C-RDC & /"MC simulation of ",i4," residues,",i6," atoms")') C-RDC & nres, natoms endif C assume that the number of residues is the number of molucules! dens = dble(nres)/boxvol call gensub(ierror) if (ierror.ne.0) return C set nstp if (nmove.ne.0) then if (nmove .gt. nres) then C-RDC write(iout,'(/" NMOVE EXCEEDS THE NUMBER OF RESIDUES:", C-RDC & /" RESET NMOVE TO THE NUMBER OF RESIDUES = ",I4)') C-RDC & nres nmove = nres endif nstp = nres/nmove i = mod(nres,nmove) if (i.gt.0) nstp = nstp + 1 elseif (nsubstep.ne.0) then nstp = min(nsub,nsubstep) nsubmv = (nsub-1)/nstp + 1 elseif (nsubmove.ne.0) then if (nsubmove.ge.nsub) then nstp = 1 nsubmv = nsub else nstp = nsub/nsubmove i = mod(nsub,nsubmove) if (i.gt.0) nstp = nstp + 1 nsubmv = nsubmove endif endif C-RDC C write some data C iensemb = 0 => NVT C iensemb = 1 => NPT if (iensemb.eq.0) then numstep = nstp*(nmc-1)+1 elseif (iensemb.eq.1) then numstep = (nstp+1)*(nmc-1)+1 endif if (wrtmc) then C-RDC write(iscr,'(/"simulate for ",i10," cycles, ",i10," steps.", C-RDC & //"*******************************", C-RDC & "********************")') C-RDC & nmc, numstep endif if (smartmc) then C-RDC write(iout, C-RDC & '(/" SMART MC SIMULATION", C-RDC & /" ------------------", C-RDC & /" NUMBER OF MC CYCLES: ",13X,I14)') nmc else C-RDC write(iout, C-RDC & '(/" MC SIMULATION",/" -------------", C-RDC & /" NUMBER OF MC CYCLES: ",13X,I14)') nmc endif if (nsubmove.eq.0) then C-RDC write(iout, C-RDC & '(" [NUMBER OF MC STEPS:",14X,I14, " ]")') C-RDC & numstep else C-RDC write(iout, C-RDC & '(" [NUMBER OF MC STEPS:",14X,I14," APPROX. ]")') C-RDC & numstep endif if (iensemb.eq.1) then C-RDC write(iout, C-RDC & '(" ENSEMBLE: ",35X,"NPT", C-RDC & /" PRESSURE: ",25X,F12.6," BAR")') C-RDC & press*1602190.0 else C-RDC write(iout, C-RDC & '(" ENSEMBLE: ",35X,"NVT", C-RDC & /" VOLUME: ",25X,F12.5," A^3", C-RDC & /" NUMBER DENSITY: ",20X,F12.5," A^-3")') C-RDC & boxvol, dens endif C-RDC write(iout, C-RDC & '(" MC-TEMPERATURE: ",23X,F9.3," K", C-RDC & /" DR: ",22X,F9.6," A", C-RDC & /" DANGLE ",22X,F9.6," DEG")') C-RDC & tempmc,drmc,dangle*degree C-RDC if (iensemb.eq.1) write(iout, C-RDC & '(" DBOX (1 DIMENSION): ",19X,F9.6," A")') C-RDC & drboxmc if (smartmc) then smartaa = (23.061*smarta)/boltzt C-RDC write(iout, C-RDC & '(" SMARTA: ",19X,F9.6," A^2")') smartaa endif C-RDC write(iout, C-RDC & '(" NUMBER OF MOVES PER MC STEP: ",I14, C-RDC & /" SEED FOR RANDOM NUMBER GENERATOR: ",I14, C-RDC & /" COORDINATES WILL BE WRITTEN EVERY",6X,I9,"-TH CYCLE", C-RDC & /" PAIRLIST WILL BE UPDATED EVERY",9X,I9,"-TH CYCLE")') C-RDC & nmove,iseed,nwrt,nupdate if (nadjst .ne. 0) then C-RDC write(iout, C-RDC & '(" DESIRED ACCEPTED/REJECTED RATIO:",8X,F8.3, C-RDC & /" ADJUST DX VALUES EVERY ",17X,I8, C-RDC & "-TH CYCLE")') ratio,nadjst endif C-RDC write(iout, C-RDC & '(/" CONVERGENCE CRITERIA:", C-RDC & /" A. WHEN COORDINATE/CHARGE FILES ARE WRITTEN", C-RDC & /" FOR ENERGY",20X,E10.3," EV", C-RDC & /" FOR DENSITY MATRIX",12X,E10.3, C-RDC & /" B. FOR ALL OTHER MC STEPS", C-RDC & /" FOR ENERGY",20X,E10.3," EV", C-RDC & /" FOR DENSITY MATRIX",12X,E10.3,/"")') C-RDC & eecrt, dencrt, eecrtmc, dencrtmc if (iensemb.eq.1) then C-RDC write(iout, C-RDC & '( //" -) MONTE CARLO ACCEPTANCE/REJECTION IS BASED ON", C-RDC & " THE ENTHALPY", C-RDC & /" (INDICATED WITH AN *).", C-RDC & //" -) AN -A- IN THE STATUS FIELD MEANS THAT THE", C-RDC & " CONFIGURATION WAS ACCEPTED,")') C-RDC write(iout,'(" -R- THAT THE CONFIGURATION WAS REJECTED,", C-RDC & /" -W- THAT THE COORDINATES/CHARGES", C-RDC & " FOR THAT STEP HAVE ", C-RDC & /" BEEN WRITTEN TO FILE,", C-RDC & /" -B- THAT THE BOX VOLUME WAS CHANGED.", C-RDC & //" -) REPORTED DENSITIES ARE NUMBER DENSITIES.")') else C-RDC write(iout, C-RDC & '( //" -) MONTE CARLO ACCEPTANCE/REJECTION IS BASED ON", C-RDC & " THE TOTAL ENERGY", C-RDC & /" (INDICATED WITH AN *).", C-RDC & //" -) AN -A- IN THE STATUS FIELD MEANS THAT THE", C-RDC & " CONFIGURATION WAS ACCEPTED,")') C-RDC write(iout,'(" -R- THAT THE CONFIGURATION WAS REJECTED,", C-RDC & /" -W- THAT THE COORDINATES/CHARGES", C-RDC & " FOR THAT STEP HAVE ", C-RDC & /" BEEN WRITTEN TO FILE.")') endif if (index(keywrd,'DIRECT').eq.0) then C-RDC write(iout,'(/" -) A DIRECT CALCULATION IS PERFORMED FOR", C-RDC & " EVERY STEP", C-RDC & /" EVEN THOUGH THIS WAS NOT SPECIFIED IN THE", C-RDC & " KEYWORDS")') endif C-RDC write(iout,'( C-RDC & /"-----------------------------------------------", C-RDC & "--------------------"//)') C adjust temperature for instanteneous pressure calculation tempmc = 1.38066D2*tempmc denstemp = dens*tempmc movedres = nres do i=1,nres ires(i) = i enddo 1000 return end CC---------------------------------------------------------------------CC subroutine mcupdate(imc,wchar,bchar,npicked,picked, & nsubpicked, subpicked, newsub, wrtstuff,nstp,ierror) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" character wchar, bchar logical newsub, wrtstuff, picked, subpicked dimension picked(maxres), subpicked(maxsub) CC------------------------------------------------------------------CC C see if data should be written to file and set convergence crit. ii = mod(imc-1,nwrt) if ((ii.eq.0).or.(imc.eq.nmc)) then wrtstuff = .true. wrtedmx = .true. else wrtstuff = .false. wrtedmx = .false. endif C set the convergence criteria to the MC-values eecrt = eecrtmc dencrt = dencrtmc wchar = ' ' bchar = ' ' C keep track of which residues are moved nscf = 0 do 76 ii=1,nres picked(ii) = .false. 76 continue do ii=1,nsub subpicked(ii) = .false. enddo npicked = 0 nsubpicked = 0 C see if new subsetting is needed ii = mod(imc-1,nupdate) if (ii .eq. 0) then newsub = .true. call gensub(ierror) if (ierror.ne.0) return if (nsubmove.ne.0) then C . calculate nstp if (nsubmove.ge.nsub) then nstp = 1 nsubmv = nsub else nstp = nsub/nsubmove i = mod(nsub,nsubmove) if (i.gt.0) nstp = nstp + 1 nsubmv = nsubmove endif elseif (nsubstep.ne.0) then nstp = min(nsub,nsubstep) nsubmv = (nsub-1)/nstp + 1 endif endif C check if adjustment to desired ratio has to be made if (nadjst .ne. 0) then if (mod(imc-1,nadjst).eq.0) then rr = dble(naccpt)/dble(ntry) if (rr .gt. ratio) then drmc = drmc*1.05 dangle = dangle*1.05 drboxmc = drboxmc*1.05 smarta = smarta*1.05 smartaa = (23.061*smarta)/boltzt smarta2 = 2.0*smartaa else drmc = drmc*0.95 dangle = dangle*0.95 drboxmc = drboxmc*0.95 smarta = smarta*0.95 smartaa = (23.061*smarta)/boltzt smarta2 = 2.0*smartaa endif C-RDC write(iout,'(/" RATIO ADJUSTMENT, NEW VALUES:", C-RDC & " (UNITS RESP A AND DEG)", C-RDC & /" DR ",F9.6, C-RDC & ", DANGLE ",F9.6, C-RDC & /" DBOX ",F9.6)') C-RDC & drmc,dangle*degree,drboxmc if (smartmc) then C-RDC write(iout,'(" SMARTA ",F9.6," A^2")') smartaa endif if (wrtmc) then C-RDC write(iscr,'(/" ratio adjustment, new values:", C-RDC & " (units resp A and deg)", C-RDC & /" dr ",F9.6, C-RDC & ", dangle ",F9.6, C-RDC & /" dbox ",F9.6)') C-RDC & drmc,dangle*degree,drboxmc if (smartmc) then C-RDC write(iscr,'(" smarta ",F9.6," A^2")') smartaa endif endif endif endif end CC----------------------------------------------------------------CC subroutine mcclose implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" character*20 name, namemcf C some closing statements C-RDC write(iout, C-RDC & '(//" MC RESTART INFORMATION: ", C-RDC & /" ----------------------", C-RDC & //" SEED = ",I18, /" DANGLE = ",F18.8, C-RDC & /" DR = ",F18.8)') C-RDC & iseed,dangle*degree,drmc if (smartmc) then smartaa = (23.061*smarta)/boltzt C-RDC write(iout,'(" SMARTA = ",F18.8)') smartaa endif C-RDC if (iensemb.eq.1) write(iout,'(" DBOX = ",F18.8, C-RDC & /" BOXDIMENSIONS ARE ",F18.8,1X,F18.8,1X,F18.8)') C-RDC & drboxmc,dbox(1),dbox(2),dbox(3) C-RDC write(iout, C-RDC & '(/" FINAL QUANTITIES:", C-RDC & /" ---------------", C-RDC & //" NSTEP = ", i10, C-RDC & /" NACCEPT = ", i10, C-RDC & /" RATIO = ", F10.3)') C-RDC & ntry, naccpt,dble(naccpt)/dble(ntry) if (nadjst .ne. 0) then C-RDC write(iout,'(" (DESIRED RATIO = ",F10.3,")")') C-RDC & ratio endif C-RDC write(iout, C-RDC & '(/" MINIMIZED ENERGIES (MC STEP ",i8,"):", C-RDC & //" ELECTRONIC ENERGY = ",F18.8," eV", C-RDC & /" CORE-CORE REPULSIONS = ",F18.8," eV")') C-RDC & nstepmc,electmc,ecoremc C-RDC if (pme) write(iout,'(" PME DIRECT ENERGY = ",F18.8," eV", C-RDC & /" PME RECIPROCAL ENERGY = ",F18.8," eV", C-RDC & /" PME SELF ENERGY = ",F18.8," eV", C-RDC & /" CLAS. COULOMB ENERGY: = ",F18.8," eV", C-RDC & /" TOTAL LONGE RANGE = ",F18.8," eV")') C-RDC & pmedirmc,pmerecmc,pmeselfmc,ecoulmc,elrmc C-RDC write(iout,'(" TOTAL ENERGY = ",F18.8," eV", C-RDC & /" HEAT OF FORMATION = ",F18.8," KCAL/MOL")') C-RDC & etotmc,eheatmc C-RDC if (iensemb.eq.1) write(iout,'(" ENTHALPY = ", C-RDC & F18.8," eV")') enthmc if (myid.eq.0) then close(iunit(7)) close(iunit(8)) if (gradient) close(iunit(10)) else C close(iunit(11)) close(iunit(17)) close(iunit(18)) if (gradient) close(iunit(20)) endif call wrtdiv(nmc) if (myid.eq.0) then name = fname(4) else name = namemcf(fname(irst),myid,ext(irst)) endif C-RDC write(iout,'(/" SAVED MINIMIZED CONFIGURATION AS DIVCON ", C-RDC & "INPUT FILE TO ",A20)') name end