! ****************************************************************
! Trial version implementing the semi-local conformational update
! BGS (Biased Gaussian Steps). This file presently contains the
! functions initlund, bgsposs and bgs.
!
! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
!                      Jan H. Meinke, Sandipan Mohanty
! ****************************************************************


! Checks if it is possible to perform a BGS update starting at the
! variable indexed ipos. Calls: none.
      logical function bgsposs(ips)
      include 'INCL.H'
      include 'incl_lund.h'
      integer jv, ips, iaa, nursvr, nnonfx, i

      logical ians

      jv=idvr(ips)
      iaa=nursvr(jv)
      ians=.true.
!      print *,'evaluating bgs possibility for ',ips,nmvr(jv)
      if (nmvr(jv).ne.'phi') then
!         print *,'bgs not possible because variable name is ',nmvr(jv)
         ians=.false.
      else if (iaa.gt.(irsml2(mlvr(jv))-3)) then
!         print *,'bgs impossible, residue too close to end'
!         print *,'iaa = ',iaa,' end = ',irsml2(mlvr(jv))
         ians=.false.
      else
         nnonfx=0
         do i=iaa,iaa+3
            if (iphi(i).gt.0) then
               if (.not.fxvr(iphi(i))) then
                  nnonfx=nnonfx+1
               endif
            endif
            if (.not.fxvr(ipsi(i))) then
               nnonfx=nnonfx+1
            endif
         enddo
         if (nnonfx.lt.6) then
!            print *,iaa,'bgs impossible because ndof = ',nnonfx
            ians=.false.
         endif
      endif
      if (ians) then
!         print *,'bgs is possible for angle ',ips,jv,nmvr(jv)
      endif
      bgsposs=ians
      return
      end

! Biased Gaussian Steps. Implements a semi-local conformational update
! which modifies the protein backbone locally in a certain range of
! amino acids. The 'down-stream' parts of the molecule outside the
! region of update get small rigid body translations and rotations.
!
! Use the update sparingly. It is rather local, and is not of great
! value if we need big changes in the conformation. It is recommended
! that this update be used to refine a structure around a low energy
! candidate structure. Even at low energies, if you always
! perform BGS, the chances of coming out of that minimum are small.
! So, there is a probability bgsprob, which decides whether BGS or the
! normal single angle update is used.
!
! Calls: energy, dummy (function provided as argument), addang, (rand)
!

      integer function bgs(eol1,dummy)
      include 'INCL.H'
      include 'incl_lund.h'
      double precision grnd, xiv, bv, ab, rv, dv, a, sum, p, r1, r2
      double precision ppsi, wfw, addang, enw, energy, wbw, rd, delta
      double precision dummy, eol1

      integer ivar, i, nph, jv, ia, nursvr, icurraa, j, k, l

      external dummy
      dimension xiv(8,3),bv(8,3),rv(3,3),dv(3,8,3)
      dimension ab(8), A(8,8),p(8),ppsi(8)
      double precision ovr(mxvr)
! Initialize
!      print *,'using BGS on angle ',nmvr(idvr(ivar))
      if (bgsnvar.eq.0) then
         bgs=0
         goto 171
      endif
      ivar=1+grnd()*bgsnvar
      do i=1,8
         iph(i)=-50000
         dph(i)=0
      enddo
      nph=0
      jv=idvr(ivar)
      ia=nursvr(jv)
! Get BGS matrices based on coordinates of atoms in 4 amino acids
      do i=1,4
         icurraa=ia+i-1
         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
            nph=nph+1
            xiv(nph,1)=xat(iCa(icurraa))
            xiv(nph,2)=yat(iCa(icurraa))
            xiv(nph,3)=zat(iCa(icurraa))
            bv(nph,1)=xiv(nph,1)-xat(iN(icurraa))
            bv(nph,2)=xiv(nph,2)-yat(iN(icurraa))
            bv(nph,3)=xiv(nph,3)-zat(iN(icurraa))
            ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
     &           +bv(nph,3)*bv(nph,3)
            iph(nph)=iphi(icurraa)
         endif
         if (.not.fxvr(ipsi(icurraa))) then
            nph=nph+1
            xiv(nph,1)=xat(iC(icurraa))
            xiv(nph,2)=yat(iC(icurraa))
            xiv(nph,3)=zat(iC(icurraa))
            bv(nph,1)=xiv(nph,1)-xat(iCa(icurraa))
            bv(nph,2)=xiv(nph,2)-yat(iCa(icurraa))
            bv(nph,3)=xiv(nph,3)-zat(iCa(icurraa))
            ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
     &           +bv(nph,3)*bv(nph,3)
            iph(nph)=ipsi(icurraa)
         endif
      enddo
      rv(1,1)=xat(iCa(ia+3))
      rv(1,2)=yat(iCa(ia+3))
      rv(1,3)=zat(iCa(ia+3))
      rv(2,1)=xat(iC(ia+3))
      rv(2,2)=yat(iC(ia+3))
      rv(2,3)=zat(iC(ia+3))
      rv(3,1)=xat(iC(ia+3)+1)
      rv(3,2)=yat(iC(ia+3)+1)
      rv(3,3)=zat(iC(ia+3)+1)

      do i=1,3
         do j=1,nph
            dv(i,j,1)=(1.0/ab(j))*(bv(j,2)*(rv(i,3)-xiv(j,3))-
     &           bv(j,3)*(rv(i,2)-xiv(j,2)))
            dv(i,j,2)=(-1.0/ab(j))*(bv(j,1)*(rv(i,3)-xiv(j,3))-
     &           bv(j,3)*(rv(i,1)-xiv(j,1)))
            dv(i,j,3)=(1.0/ab(j))*(bv(j,1)*(rv(i,2)-xiv(j,2))-
     &           bv(j,2)*(rv(i,1)-xiv(j,1)))
         enddo
      enddo
      do i=1,nph
         do j=i,nph
            A(i,j)=0
            do k=1,3
               do l=1,3
                  A(i,j)=A(i,j)+dv(k,i,l)*(dv(k,j,l))
               enddo
            enddo
            A(i,j)=bbgs*A(i,j)
            if (i.eq.j) then
               A(i,j)=A(i,j)+1
            endif
            A(i,j)=0.5*abgs*A(i,j)
         enddo
      enddo
      do i=1,nph
         do j=i,nph
            sum=A(i,j)
            do k=i-1,1,-1
               sum=sum-A(i,k)*A(j,k)
            enddo
            if (i.eq.j) then
               p(i)=sqrt(sum)
            else
               A(j,i)=sum/p(i)
            endif
         enddo
      enddo
! Generate 8 Gaussian distributed small random angles
      do i=1,8,2
         r1=grnd()
!        In the rare event that this random number is 0, just take the next
         if (r1.le.0) r1=grnd()
         r1=sqrt(-log(r1))
         r2=grnd()
         ppsi(i)=r1*cos(pi2*r2)
         ppsi(i+1)=r1*sin(pi2*r2)
      enddo
      do i=1,nph
         dph(i)=0
      enddo
! Solve lower triangular matrix to get dphi proposals
      do i=nph,1,-1
         sum=ppsi(i)
         do k=i+1,nph
            sum=sum-A(k,i)*dph(k)
         enddo
         dph(i)=sum/p(i)
      enddo
! Calculate intrinsic (non-Boltzmann) weight for forward process
!      print *,'calculating intrinsic weight for forward process'
      sum=0
      do i=1,nph
         sum=sum+ppsi(i)*ppsi(i)
      enddo
      wfw=exp(-sum)
      do i=1,nph
         wfw=wfw*p(i)
      enddo

! Reconstruct chain and calculate new energy
!      print *,'going to assign changes to the chain'
      ovr = vlvr
      do i=1,nph
         vlvr(iph(i))=addang(vlvr(iph(i)),dph(i))
      enddo
      enw = energy()
! Calculate weight for reverse process for detail balance
!      print *,'proceeding to calculate weight for the reverse process'
      nph=0
      do i=1,4
         icurraa=ia+i-1
         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
            nph=nph+1
            xiv(nph,1)=xat(iCa(icurraa))
            xiv(nph,2)=yat(iCa(icurraa))
            xiv(nph,3)=zat(iCa(icurraa))
            bv(nph,1)=xiv(nph,1)-xat(iN(icurraa))
            bv(nph,2)=xiv(nph,2)-yat(iN(icurraa))
            bv(nph,3)=xiv(nph,3)-zat(iN(icurraa))
            ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
     &           +bv(nph,3)*bv(nph,3)
            iph(nph)=iphi(icurraa)
         endif
         if (.not.fxvr(ipsi(icurraa))) then
            nph=nph+1
            xiv(nph,1)=xat(iC(icurraa))
            xiv(nph,2)=yat(iC(icurraa))
            xiv(nph,3)=zat(iC(icurraa))
            bv(nph,1)=xiv(nph,1)-xat(iCa(icurraa))
            bv(nph,2)=xiv(nph,2)-yat(iCa(icurraa))
            bv(nph,3)=xiv(nph,3)-zat(iCa(icurraa))
            ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
     &           +bv(nph,3)*bv(nph,3)
            iph(nph)=ipsi(icurraa)
         endif
      enddo
      rv(1,1)=xat(iCa(ia+3))
      rv(1,2)=yat(iCa(ia+3))
      rv(1,3)=zat(iCa(ia+3))
      rv(2,1)=xat(iC(ia+3))
      rv(2,2)=yat(iC(ia+3))
      rv(2,3)=zat(iC(ia+3))
      rv(3,1)=xat(iC(ia+3)+1)
      rv(3,2)=yat(iC(ia+3)+1)
      rv(3,3)=zat(iC(ia+3)+1)

      do i=1,3
         do j=1,nph
            dv(i,j,1)=(1.0/ab(j))*(bv(j,2)*(rv(i,3)-xiv(j,3))-
     &           bv(j,3)*(rv(i,2)-xiv(j,2)))
            dv(i,j,2)=(-1.0/ab(j))*(bv(j,1)*(rv(i,3)-xiv(j,3))-
     &           bv(j,3)*(rv(i,1)-xiv(j,1)))
            dv(i,j,3)=(1.0/ab(j))*(bv(j,1)*(rv(i,2)-xiv(j,2))-
     &           bv(j,2)*(rv(i,1)-xiv(j,1)))
         enddo
      enddo
      do i=1,nph
         do j=i,nph
            A(i,j)=0
            do k=1,3
               do l=1,3
                  A(i,j)=A(i,j)+dv(k,i,l)*(dv(k,j,l))
               enddo
            enddo
            A(i,j)=bbgs*A(i,j)
            if (i.eq.j) then
               A(i,j)=A(i,j)+1
            endif
            A(i,j)=0.5*abgs*A(i,j)
         enddo
      enddo
      do i=1,nph
         do j=i,nph
            sum=A(i,j)
            do k=i-1,1,-1
               sum=sum-A(i,k)*A(j,k)
            enddo
            if (i.eq.j) then
               p(i)=sqrt(sum)
            else
               A(j,i)=sum/p(i)
            endif
         enddo
      enddo
      do i=1,nph
         ppsi(i)=p(i)*dph(i)
         do j=i+1,nph
            ppsi(i)=ppsi(i)+A(j,i)*dph(j)
         enddo
      enddo
      sum=0
      do i=1,nph
         sum=sum+ppsi(i)*ppsi(i)
      enddo
      wbw=exp(-sum)
      do i=1,nph
         wbw=wbw*p(i)
      enddo
!     Acceptance condition (includes additional weight)
      rd=grnd()
!      print *,'generated  selection random number ',rd
!      print *,'wfw/wbw = ',wfw/wbw
      rd=-log(rd*wfw/wbw)
!      print *,'modified rd = ',rd
!      print *,'before calculating energy change'
      delta = dummy(enw)-dummy(eol1)
!      print *,'delta = ',delta
!       call outpdb(0,'after.pdb')
!      print *,'after outpdb for after.pdb'
!      do i=1,nph
!         print *,'BGS>',i,iph(i),vlvr(iph(i)),dph(i)
!      enddo
      if (rd.ge.delta) then
!     accept
         eol1=enw
         bgs=1
!         print *,'BGS move accepted'
      else
!     reject
         vlvr = ovr
!          enw=energy()
!          if (abs(enw-eol1).gt.0.000001) then
!             write (logString, *) 'rejected bgs move: energy change :',eol1,enw
!          endif
!         write (logString, *) 'rejected bgs move: ',eol1,enw,wfw,wbw,rd
         bgs=0
      endif
 171  continue
      return
      end