program hybrid_exp3
c--------------------------------------------------------------------
c     program hybrid and variants, by Michael K. Gilson
c     Copyright 1992, all rights reserved.
c     For non-commercial use only.  
c     Academic distribution permitted, but please notify me upon
c     distribution.  Also, please inform me of any bugs you may find,
c     and of any significant enhancements.  Thank you.
c     This code comes with no guarantee.
c---------------------------------------------------------------------
c Hybrid of tanford and roxby mean field method and exact method 
c for titration curves, using combined reduced site and cluster methods.  
c Reference: Gilson,M.K. Multiple-Site Titration and Molecule Modelling:
c Two Rapid Methods for Computing Energies and Forces for Ionizable
c Groups in Proteins.  Proteins: Structure, Function and Genetics, 
c 15:266-282, 1993.
c 
c This version includes relaxation to avoid oscillation in the iterative 
c loop.  That is, it uses the average of the last two iterations in theta 
c to compute effect of another cluster on current cluster, in subroutine 
c thetacluster.  (I have fixed the average to a 50-50 weighting, but in 
c some cases another weighting might be preferable.)
c
c This version loops over phs, but this loop should be easy to
c remove, if desired.  
c
c Each "fixed" group is treated  as a separate cluster, thus satisfying 
c Bashford and Karplus's idea of using the mean field approximation for 
c the fixed groups.  Automatically adjusts energy cutoff for cluster 
c separation, to force maximum cluster size to be observed.  If
c one cluster is too big for a given Gc, Gc is increased just for
c this cluster, until its size is within the limit.  However, other
c clusters may fall within the allowed size with smaller Gc values.
c 
c Output includes average of purely electrostatic energy (Gelec).
c
c This code at one time offered the option of looping over
c ionizable groups in the forward order or backward order, but
c I have commented out the option... now it only goes forward.
c I have also commented out some timing printout (maximum time per
c calculation at a given pH, and average time for all the pHs).
c
c Some options and printouts have been commented out, but not deleted 
c to preserve them as options which can be revived.
c
c This code uses a recursive subroutine call (subroutine tree) which will not
c work on all systems. 
c
c A reported pKa of 999. is not an error; it is a flag indicating
c that the range of pHs you requested did not bracket the pKa of the group.
c
c A "Largest cluster size" of 999 also is not an error; it indicates
c that all groups were "fixed".  This usually occurs at very high
c or low pHs.
c-------------------------------------------------------------
c working version as of 12/16/92
c hybrid of tanford and roxby mean field method and 
c exact method for titration curves. uses reduced site and cluster
c methods
c version includes  relaxation to avoid
c oscillation.  uses  average oflast two theta values in
c getting effect of another cluster on current cluster,
c in subroutine thetacluster.  (I have fixed the average to
c a 50-50 average, but in some cases another weighting could be 
c prefereable.  
c  loops over phs, does timing
c this version treats each fixed group as a separate cluster,
c thus satisfying bashford's idea of doing mean field for
c the fixed groups
c current version automatically adjusts energy cutoff for
c cluster separation, based upon some maximum cluster size.
c also...returns average of purely electrostatic energy
c
c does not recompute all cluster sizes if one is too big: uses
c new cluster subroutine which keeps restarting cluster of a certain
c residue until it is small enough.  
c for some reason, this code offered the option of looping over
c ionizable groups in the forward order or backward order, but
c I have commented out the option... now it only goes forward.
c I have also commented out some timing printout (max time per
c calc at a given pH, and average time for all the pHs).

      implicit double precision (a-h,o-z)
      parameter (factor=1./(2.303*.6))
      parameter (beta=1./.6)

      logical fixed(500)
      dimension pka0(500),pkai(500),pka(500),phs(500)
      dimension thetaold(500),thetanew(500),thetaolder(500)
      dimension thetastore(500,500),tstore(500)
      dimension g(500,500),gbig(500), pkapp(500)
      integer clust(500,500),nset(500),z(500),maxc(2),ibig(500)
      integer pos, neg
      real*4 t(2),dtime
      character*80 infile,outfile
      data maxc(1) /500/
      tmax = -99
      tav = 0.
      nrun = 0
c      print *, 'delta theta for convergence, weight factor:'
c      accept *, conv,wt
      conv = .000001
      wt = 0.5
      print *, 'preferred max cluster size (suggested 10):'
      read (5,*) maxc(2)
      nmax=2
c preprocessing for clusters larger than some maximum (here 10)
      if (maxc(2).le.maxc(1)) then
         maxc(1) = maxc(2)
         nmax = 1
      end if

      print *, 'fixed cutoff theta (suggested 0.05):'
      read (5,*) tcut

      print *, 'phstart,phstop,phstep:'
      read (5,*) ph1,ph2,phstep
      nph = (ph2-ph1)/phstep +1

c      print *, 'run groups forward or backward : 0 or 1'
c      accept *, idirection

 50   continue
c read in data file
      print *, 'input file:'
      read (5,991) infile
 991  format (a80)
      print *, 'outfile:'
      read (5,991) outfile
      open (unit=10,file=infile)
      open (unit=11,file=outfile)
      read (10,*) ngroup
c      print *, 'ngroup:',ngroup
      print *, 'ignore self energies? 1 for yes (not suggested):'
      read (5,*) iself
      do i = 1, ngroup
         read (10,*) pka0(i),z(i),g(i,i)
         if (z(i).lt.0) neg = neg+z(i)
         if (z(i).gt.0) pos = pos+z(i)
         if (iself.ne.1) then
            pkai(i) = pka0(i) -factor*z(i)*g(i,i)
         else
            pkai(i) = pka0(i)
         end if
c            print *,  'pkai:',pkai(i)
         do j = i+1,ngroup
            read (10,*) g(i,j)
            g(j,i) = g(i,j)
         end do
      end do
      close (10)
      print *, 'Greatest possible neg charge:', neg
      print * ,'Greatest possible pos charge:', pos
      idirection=0
      if (idirection.eq.0) then
         i1 = 1
         i2 = ngroup
         i3  = 1
      else
         i1 = ngroup
         i2 = 1
         i3 = -1
      end if
      write (6,9998)
 9998 format 
     &(1x, '   ph   N clusts   Lrgst clust    Gc(max)')

      write (11,9993)nph
 9993 format (1x,i6)

      write (11,9991) 
 9991 format (1x,
     &'  ph     mean q     Gion       dGion      ',
     &'CPU(s)   gcut   tcut    Gelec')

      do iph = 1,nph
         nrun = nrun+1
c loop over states, summing contributions to perturbed (numerator)
c and unperturbed (denominator) ionization energies.
         ph = ph1+(iph-1)*phstep
         phs(iph) = ph
         time= dtime(t)

c assign initial charge states    
c start with everything ionized in this version  
c thetanew, thetaold = thetanew
         
         do i = i1, i2,i3
c            call charger(i,pkai,ph,z,thetaold)
            thetaold(i) = 1.
            thetanew(i) = thetaold(i)
            thetaolder(i) = thetaold(i)
         end do

c assign fixed charges for this ph
         call fixer(nfixed,fixed,thetaolder,thetaold,
     &        thetanew, pkai,z,ph,g,ngroup,tcut)

         do imaxcluster = 1,nmax
            maxcluster = maxc(imaxcluster)

c assign clusters for this ph (depends on fixed groups)
 200     continue
         call clust2(ngroup,fixed,clust,nset, nc,g,maxcluster,
     &        gcmin,gcmax)
         nsmax = -999

         do i  =1, nc
            nsmax = max(nset(i),nsmax)
c            print *, 'nset:', nset(i)
c            do ii = 1,nset(i)
c               print *, clust(ii,i)
c            end do
         end do

         if (imaxcluster.eq.nmax) then
            write (6,9999) ph, nc, nsmax, gcmax
 9999       format (1x, f7.3, 2x, i4, 7x, i4,10x, f8.5)
         end if
c stick fixed residues in individual clusters:
         do if = 1, ngroup
            if (fixed(if)) then
               nc = nc+1
               nset(nc) = 1
               clust(1,nc) = if
            end if
         end do

       nit = 1
 100   continue
c loop over clusters
       gionsum = 0.
       gelecsum = 0.
       do icluster = 1, nc
         
c do exact computation of charges states within this cluster, using
c mean field approx for influence of fixed charges and of charges
c outside the cluster. use thetaold arrays in calc, generate thetanew
c arrays

c NOTE THAT thetacluster's gion overestimates inter-cluster
c electrostatic interactions by a factor of 2 relative to the
c intended approximation; this problem is compensated for subsequently:

          call thetacluster(icluster,clust,nset,thetaolder,thetaold,
     &         thetanew,g,ph,pkai,z,ngroup,fixed,nc,gion,wt,gelec)

          gelecsum = gelecsum + gelec
          gionsum = gionsum + gion
       end do

c check for convergence (thetanew vs thetaold)
c also allow convergence on first pass if there's only
c one cluster; so no iterative correction required:
       d1max = -1
       d2max = -1
       do igroup = 1, ngroup
          d1 = abs(thetanew(igroup)-thetaold(igroup))
          d2 = abs(thetaold(igroup)-thetaolder(igroup))
          if (d1.gt.d1max) then
             d1max = d1
             i1max = igroup
          end if
          if (d2.gt.d2max) then
             d2max = d2
             i2max = igroup
          end if

          if ((d1.gt.conv )
     &         .and. nc.ne.1) then
c     if continuing, set thetaold = thetanew, and goto 100
c             print *, 'igroup,new',
c     &            igroup,thetanew(igroup)
             do jgroup = 1, ngroup
                thetaolder(jgroup) = thetaold(jgroup)
                thetaold(jgroup) = thetanew(jgroup)
             end do
             nit = nit+1
c             if (nit.gt.20) print *, 'not converging'
             goto 100
          end if
       end do
c convergence achieved for this cluster size
c do next cluster size
       end do
c finally converged
       
c get ionization energy exactly for unperturbed pkas:
       tstore (iph) = dtime(t)
       tmax = max(tmax,tstore(iph))
       tav = tav+tstore(iph)

       call gexact(pka0,ph,z,ngroup,gex0)
c compute mean interaction among clusters and subtract from
c gionsum to correct for factor of two error in this term
c incurred by thetacluster:
c also look for groups with inter-cluster interactions gt 1kcal/mole, say.
       gsum = 0.
       nbig = 0
       do ic =1,nc
          do iset = 1, nset(ic)
             ig = clust(iset,ic)
c     careful to single-count interactions
             do jc = ic+1, nc
                do jset = 1,nset(jc)
                   jg = clust(jset,jc)
                   if (g(ig,jg).gt.1.5 .and. 
     &                  .not.fixed(ig) .and..not.fixed(jg)) then
                      nbig = nbig+1
                      ibig(nbig) = ig
                      gbig(nbig) = g(ig,jg)
                   end if
                   gsum = gsum + g(ig,jg)*thetanew(ig)*
     &                  thetanew(jg)*z(ig)*z(jg)
                end do
             end do
          end do
       end do
c       print *, 'possible problem groups for ph:', ph
c       do i= 1, nbig
c          print *, ibig(i),gbig(i)
c       end do

c     print *, 'Cluster-cluster interactions: ',gsum
       gelecsum = gelecsum + gsum
       qsum = 0.
c compute purely electrostatic mean self energy of charges,too.
       gself =0
       do i = 1, ngroup
          qsum = qsum + thetanew(i)*z(i)
          thetastore(i,iph) = thetanew(i)
          gself = gself + thetanew(i)*g(i,i)
       end do
       gelecsum = gelecsum + gself
c write info:
       if (gcmax.lt.0) gcmax = 9.999
         write (11,9992) ph, qsum, gionsum-gsum, gionsum-gsum-gex0,
     &      tstore(iph), gcmax,tcut,gelecsum

 9992    format (1x, f7.2, f8.2,f12.2,f10.2,f12.6,1x,f6.3,1x,f6.3,
     &        f10.2)

c end ph loop
      end do

c write charge info      

      write (11,9996) ngroup
 9996 format (1x,i6)

c estimate and write out apparent pkas.
      do igroup =1,ngroup
         pkapp(igroup) = 999.
      end do

      do iph = 1,nph - 1
         do i=1,ngroup
            if ( (thetastore(i,iph) .lt. 0.500 .and. 
     &           thetastore(i,iph+1) .gt. 0.500) .or.
     &           (thetastore(i,iph) .gt. 0.500 .and. 
     &           thetastore(i,iph+1) .lt. 0.500) ) then

               slope = (phs(iph+1)-phs(iph))/
     &              (thetastore(i,iph+1)-thetastore(i,iph))

               pkapp(i) = phs(iph) + slope*(0.500 - thetastore(i,iph))
            end if
         end do
      end do
      write (11,999) 
 999  format (1x, 'Group    pk(model)  pK(app)     pK(app) - pK(model)')
      do i=1,ngroup
         write (11,9997) i, pKa0(i),   pkapp(i), pkapp(i) - pka0(i)
      end do
 9997 format (i4, 3x, f8.3, 5x, f8.3, 8x, f8.3)

      do iph = 1,nph
         write (11,9994) phs(iph)
         do i=1,ngroup
            write (11,9995) thetastore(i,iph)*z(i)
         end do
      end do
 9994 format (1x,f7.2)
 9995 format (1x,f8.3)

      close (11)
c      print *, 'tmax:', tmax
c      print *, 'tav to date:', tav/nrun
      goto 50
      end
c---------------------------------------------------
      subroutine charger(i,pka,ph,z,theta)
      implicit double precision (a-h,o-z)
      dimension pka(500),theta(500)
      integer z(500)

      d = 10.**(z(i)*(pka(i)-ph))
      theta(i) = d/(1+d)
      return
      end
c-------------------------------------------------------
      subroutine fixer(nfixed,fixed,thetaolder,thetaold,
     &     thetanew,pkai,z,ph,g,ngroup,tcut)
      implicit double precision (a-h,o-z)
      parameter (factor=1./(2.303*.6))      
      dimension thetaold(500),pkai(500),g(500,500)
      dimension pka1(500),pka2(500)
      dimension theta1(500),theta2(500),thetanew(500)
      dimension thetaolder(500)
      logical fixed(500)
      integer z(500)

      nfixed = 0
      tctop = 1.-tcut
      do i = 1, ngroup
         fixed(i) = .false.

c     compute new pka with protein maximally positive:
         g1 =0.
         do j= 1, ngroup
c     last part of next line is trick to include only basic
c     residues, without using conditional:
            if (j.ne.i) g1= g1 + z(i)*z(j)*g(i,j) * (z(j)+1.) * .5
         end do
         pka1(i) = pkai(i) - factor*z(i)*g1
c     get theta for this case:
         call charger(i,pka1,ph,z,theta1)

c     compute new pka with protein maximally negative:
         g2 = 0.
         do j = 1, ngroup
            if (j.ne.i) g2 = g2+ z(i)*z(j)*g(i,j) * (z(j)-1.)*(-.5)
         end do

         pka2(i) = pkai(i) -factor*z(i)*g2
c     get theta for this case:         
         call charger (i,pka2,ph,z,theta2)

         if ((theta2(i).lt.tcut .and. theta1(i).lt.tcut)
     &        .or.(theta2(i).gt.tctop.and.theta1(i).gt.tctop)) then
            nfixed = nfixed +1
            fixed(i) = .true.


            thetaold(i) = theta1(i)
            thetanew(i) = thetaold(i)
            thetaolder(i) = thetaold(i)
         end if
      end do

      return
      end
c---------------------------------------------------
      subroutine clust2(ngroup,fixed,clust,nset, nc,g,maxcluster,
     &     gcmin,gcmax)
      implicit double precision (a-h,o-z)
      dimension g(500,500)
      logical used(500),fixed(500)
      integer clust(500,500),nset(500)

      do igroup = 1, ngroup
         used(igroup) = .false.
      end do
      gcmin = 1000
      gcmax = -1000

      nc = 0
      do igroup = 1, ngroup
c if group not yet in cluster:
         if (.not.used(igroup) .and. .not. fixed(igroup)) then
c start new cluster:
            nc = nc+1
            gc = 0.
 100        iset= 1
            used(igroup) = .true.
            clust(iset,nc) = igroup
c begin recursive search from this group:
            call tree(iset,nc,used,clust,g,gc,igroup,
     &           igroup+1,ngroup,fixed)
c number of groups in this cluster:
            nset(nc) = iset
            gcmin = min(gcmin,gc)
            gcmax = max(gcmax,gc)

c if cluster is too big, increase energy cutoff and split it smaller:
            if (nset(nc).gt.maxcluster) then

c reset members to "unused" status
               do iset = 1, nset(nc)
                  used(clust(iset,nc))=.false.
               end do
c increase cutoff energy
               gc = gc + .1
c recalculate this cluster:
               goto 100
            end if
         end if
      end do

      return
      end
c------------------------------------------
      subroutine tree(iset,nc,used,clust,g,gc,igroup,
     &     istart,ngroup,fixed)

      implicit double precision (a-h,o-z)
      dimension g(500,500)
      integer clust(500,500)
      logical used(500),fixed(500)

c start search with group which initiated this cluster
c this will shorten search towards the end.  
      do jgroup = istart,ngroup
         if (.not.used(jgroup) .and. .not. fixed(jgroup)) then
            if (g(igroup,jgroup).gt.gc)then
               used(jgroup) = .true.
               iset = iset+1
               clust(iset,nc) = jgroup
               call tree(iset,nc,used,clust,g,gc,jgroup,
     &              istart,ngroup,fixed)

            end if
         end if
      end do
      return
      end

c----------------------------------------------------------
      subroutine thetacluster(icluster,clust,nset,thetaolder,
     &     thetaold,thetanew,g,ph,pkai,z,ngroup,fixed,nc,gion,wt,
     &     gelec)
      implicit double precision (a-h,o-z)
      parameter (factor=1./(2.303*.6))
      parameter (beta=1./0.6)
      dimension pkai(500)
      dimension thetaold(500),thetanew(500),thetaolder(500)
      dimension g(500,500)
      integer clust(500,500),nset(500),z(500)
      integer zalph(500)
      double precision ka,kpart(500),kprod
      logical fixed(500)

      h=10.**-pH
c assign equilibrium constants using intrinsic pkas modified by
c average charges of groups in other clusters
      do iset=1,nset(icluster)
         gsumc = 0.
c index of this group:
         ig = clust(iset,icluster)
c other clusters' influence:         
         do jcluster = 1, nc
            if (jcluster.ne.icluster) then
               do jset = 1,nset(jcluster)
c index of group j:
                  jg = clust(jset,jcluster) 
c weight influence of group jg by its fractional charge
c (mean field approximation) (use average from last two iterations
c to damp possible oscillations):
                  gsumc = gsumc + g(ig,jg)*z(ig)*z(jg)
     &                 *(wt*thetaold(jg)+(1-wt)*thetaolder(jg))
               end do
            end if
         end do

         pka = pkai(ig) - factor*z(ig)*gsumc
         ka = 10.**(-pka)
         kpart(ig) = ka**(-z(ig)) * h**z(ig)
      end do
c now we have a current set of pkas, translated into kas.

      nstates = 2.**nset(icluster)

c initialize mean charges to zero
      do iset = 1,nset(icluster)
         ig = clust(iset,icluster)
         thetanew(ig) = 0.
      end do

c loop over ionization states of this cluster:
      sigma = 0.
      gelec = 0.
      do m = 1, nstates
         call zfill(zalph,nset(icluster),m)
         kprod = 1.
         gsum = 0.
c loop over groups for this state:
         do iset = 1, nset(icluster)
            ig = clust(iset,icluster)
c accumulate occupancy factor for this state:
            if (zalph(iset).eq.1) then
c include mass action effect for each group ionized in this state:
               kprod = kprod* kpart(ig)
c work on total elect energy for state m; allow different charges q=+/-1
c including j = i gives self energy term. keeping j .ge.i prevents
c double counting
c includ self energy to get correct net electrostatic energy:
               do jset = iset+1,nset(icluster)
                  jg = clust(jset,icluster)
                  gsum = gsum + zalph(jset)*g(ig,jg)*z(ig)*z(jg)
               end do
            endif
         end do
c electrostatic boltzmann factor for this state            
         bfact = exp(-gsum*beta)
c mass action weighted by boltzmann factor:
         term = bfact*kprod
c average electrostatic energy within this cluster
         gelec = gelec + term*gsum
c partition function:
         sigma = sigma + term
c accumulate ionization fractions of groups:
         do iset = 1, nset(icluster)
            if (zalph(iset) .eq.1) then
               ig = clust(iset,icluster)
c if group ionized in this state, weight charge by occupancy of state:
               thetanew(ig) = thetanew(ig)+term
            end if
         end do
c goto next state:
      end do

c get fractional charges:
      do iset = 1,nset(icluster)
         ig = clust(iset,icluster)
         thetanew(ig) = thetanew(ig)/sigma
c store charges for unperturbed run:
      end do
      gion = (-1./beta)*log(sigma)
      gelec = gelec/sigma

      return
      end
c----------------------------------------------
c fills z with all possible charge states.
      subroutine zfill(z,n,m)
      integer z(500)

      a = m-1
      do i = n-1,0,-1
c         ff = 2.**i
         ff = ishft(1,i)
c         z(i+1) = a/ff
         ia = a
         z(i+1) = ishft(ia,-i)
         a = a- ff*z(i+1)
      end do

      return
      end
         
c-------------------------------------------------------
      subroutine gexact(pk,ph,z,ngroup,gion)
      implicit double precision (a-h,o-z)
      dimension pk(500)
      integer z(500)
      parameter (RT=.6)
      gion =0.
      do i = 1, ngroup
         gion = gion - RT*log(1.+10.**(z(i)*(pk(i)-ph)))
      end do
      return
      end