! **************************************************************
!
! This file contains the subroutines: opereg,gdtgbl,gdtreg
!
! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
!                      Shura Hayryan, Chin-Ku 
! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
!                      Jan H. Meinke, Sandipan Mohanty
!
! **************************************************************

      subroutine opereg(nml)

! .......................................................................
! PURPOSE: calculate regul. energy & it's partial derivatives
!          for molecule 'nml' vs. variables 'iv'
!
!  NB: if the unit axis for an internal variable coincides with a
!      global axis (i.e. for torsion or bond length variation round
!      or along 'xrfax', respectively, and bd. angle var. round
!      'zrfax'): VdW & 14 interaction partners of moving set atoms
!      should be used for calculation, instead of the mov. sets,
!      with opposite sign.
!
!      Example: By the the way the molecule-fixed system is set up,
!               changes in Phi_1 affect atomic positions BEFORE the
!               N-C^alpha bond relatively to the space-fixed system,
!               not the moving set of Phi_1.
!
! CALLS:    gdtgbl, gdtreg
! ......................................................................

      include 'INCL.H'
      include 'INCP.H'

      dimension xfat(mxat),yfat(mxat),zfat(mxat),
     &          xfrat(mxat),yfrat(mxat),zfrat(mxat),

     &          xfvr(mxvr),yfvr(mxvr),zfvr(mxvr),
     &          xfrvr(mxvr),yfrvr(mxvr),zfrvr(mxvr)

      logical   lnb


      ntlvr=nvrml(nml)
      if (ntlvr.eq.0) then
        write (*,'(a,i4)')
     &           ' opereg> No variables defined in molecule #',nml
        return
      endif

      ix2=ixrfpt(2,nml)    ! as indicator for situation noted above

      ifivr=ivrml1(nml)         ! 1st var. &
      ilavr=ifivr+ntlvr-1       ! last var. of 'nml'

! --------------------------- initializations
      do i=ifivr,ilavr
        gdeyrg(i)=0.d0
        xfvr(i)=0.d0
        yfvr(i)=0.d0
        zfvr(i)=0.d0
        xfrvr(i)=0.d0
        yfrvr(i)=0.d0
        zfrvr(i)=0.d0
      enddo

      ii=(nml-1)*6

      do i=ii+1,ii+6
        gdeygb(i) = 0.d0
      enddo

      x1=rfpt(1,nml)     ! r_1
      y1=rfpt(2,nml)
      z1=rfpt(3,nml)

      a= gbpr(4,nml)     ! alpha
      sa = sin(a)
      ca = cos(a)

      xk = yrfax(1,nml)  ! axis K
      yk = yrfax(2,nml)
      zk = yrfax(3,nml)

      eyrg = 0.d0

      do i=iatrs1(irsml1(nml)),iatrs2(irsml2(nml))

        j = ixatp(i)
        if (j.gt.0) then

          xi = xat(i) ! position of atom in internal model
          yi = yat(i)
          zi = zat(i)

          xji = xatp(j) - xi    ! x distance between internal model and PDB
          yji = yatp(j) - yi
          zji = zatp(j) - zi

          eyrg = eyrg + xji**2 + yji**2 + zji**2  ! The regularization energy is just
                                                  ! the sum over the atom distances 
                                                  ! squared.

          dx = 2.d0 * xji   ! f = - dE/dR_i
          dy = 2.d0 * yji   ! The factor of 2 comes from the derivative
          dz = 2.d0 * zji

! =============================================== global pars.

          gdeygb(ii+1) = gdeygb(ii+1) - dx   ! d(E_ij) / d(x_i)
          gdeygb(ii+2) = gdeygb(ii+2) - dy   ! d(E_ij) / d(y_i)
          gdeygb(ii+3) = gdeygb(ii+3) - dz   ! d(E_ij) / d(z_i)

! -------------------------- r = r_i - r_1
          x = xi - x1
          y = yi - y1
          z = zi - z1

          gdeygb(ii+4) = gdeygb(ii+4) +dx*y-dy*x                      !  d(E_ij) / d(a)

          gdeygb(ii+5) = gdeygb(ii+5) +z*(dy*ca-dx*sa)+dz*(x*sa-y*ca) !  d(E_ij) / d(b)

          gdeygb(ii+6) = gdeygb(ii+6) +dx*(zk*y-yk*z)+dy*(xk*z-zk*x)  !  d(E_ij) / d(g)
     &                                +dz*(yk*x-xk*y)

! =============================================== for internal vars.

          xfat(i) = dx
          yfat(i) = dy
          zfat(i) = dz

          xfrat(i) = dy*zi-dz*yi  ! g = f x r
          yfrat(i) = dz*xi-dx*zi  !
          zfrat(i) = dx*yi-dy*xi  !

        else
          xfat(i) = 0.d0
          yfat(i) = 0.d0
          zfat(i) = 0.d0
          xfrat(i) = 0.d0
          yfrat(i) = 0.d0
          zfrat(i) = 0.d0
        endif

      enddo ! atoms

      if (tesgrd) call gdtgbl(nml)

      i1s=imsml1(nml)+nmsml(nml)  ! last mov. set of 'nml' + 1
      i1a=iadml1(nml)+nadml(nml)  ! last added var. of 'nml' + 1

      do io=ilavr,ifivr,-1  ! ______ loop over vars in desc. order

        lnb = .false.  ! = true, if situation noted above takes place

        iv=iorvr(io)       ! index,
        it=ityvr(iv)       ! type,
        ia=iatvr(iv)       ! primary mov. atom,
        ib=iowat(ia)       ! "base" of current var.

        xb=xat(ib)
        yb=yat(ib)
        zb=zat(ib)

! ---------------------------------------- axis for var.

        if (it.eq.3) then      ! torsion

          ex= xtoat(ib)
          ey= ytoat(ib)
          ez= ztoat(ib)

          if (ib.eq.ix2) lnb = .true.

        elseif (it.eq.2) then  ! b.angle

          ex= xbaat(ia)
          ey= ybaat(ia)
          ez= zbaat(ia)

          if (ib.eq.ix2) lnb = .true.

        elseif (it.eq.1) then  ! b.length

          ex=xtoat(ia)
          ey=ytoat(ia)
          ez=ztoat(ia)

          if (ia.eq.ix2) lnb = .true.

        endif

        xfiv=0.0
        yfiv=0.0
        zfiv=0.0
        xfriv=0.0
        yfriv=0.0
        zfriv=0.0

        if (.not.lnb) then

          i2s=i1s-1      ! last m.s  &
          i1s=imsvr1(iv) ! 1st m.s for var. index 'iv'

          do ims=i1s,i2s  ! __ loop over moving sets

            i1=latms1(ims)   ! 1st &
            i2=latms2(ims)   ! last mov. atom in mov. set 'ims'

            do i=i1,i2  ! __ loop over atoms i ===================

              xfiv = xfiv + xfat(i)   ! f
              yfiv = yfiv + yfat(i)   !
              zfiv = zfiv + zfat(i)   !

              xfriv = xfriv + xfrat(i)  ! g
              yfriv = yfriv + yfrat(i)  !
              zfriv = zfriv + zfrat(i)  !

            enddo  ! ... atoms i
          enddo  ! ... m.s.

          i2a=i1a-1       ! last &
          i1a=iadvr1(iv)  ! 1st 'added' var. for 'iv'

          do iad=i1a,i2a  ! loop over add. var.

            lad=ladvr(iad)

            xfiv = xfiv + xfvr(lad)
            yfiv = yfiv + yfvr(lad)
            zfiv = zfiv + zfvr(lad)
            xfriv = xfriv + xfrvr(lad)
            yfriv = yfriv + yfrvr(lad)
            zfriv = zfriv + zfrvr(lad)

          enddo

        else

          do ivw=ivwat1(ia),ivwat2(ia)  ! vdW-domains of 'ia'
            do j=lvwat1(ivw),lvwat2(ivw)  ! .. their atoms

              xfiv = xfiv - xfat(j)
              yfiv = yfiv - yfat(j)
              zfiv = zfiv - zfat(j)

              xfriv = xfriv - xfrat(j)
              yfriv = yfriv - yfrat(j)
              zfriv = zfriv - zfrat(j)

            enddo
          enddo

          do i14=i14at1(ia),i14at2(ia)  ! 1-4 partn. of 'ia'
            j=l14at(i14)

            xfiv = xfiv - xfat(j)
            yfiv = yfiv - yfat(j)
            zfiv = zfiv - zfat(j)

            xfriv = xfriv - xfrat(j)
            yfriv = yfriv - yfrat(j)
            zfriv = zfriv - zfrat(j)

          enddo

        endif

        xfvr(iv) = xfiv
        yfvr(iv) = yfiv
        zfvr(iv) = zfiv
        xfrvr(iv) = xfriv
        yfrvr(iv) = yfriv
        zfrvr(iv) = zfriv

        if (it.eq.3.or.it.eq.2) then  ! torsion,b.angle

          gdeyrg(iv)= (ey*zb-ez*yb)*xfiv+(ez*xb-ex*zb)*yfiv+
     &                (ex*yb-ey*xb)*zfiv
     &               +ex*xfriv+ey*yfriv+ez*zfriv

        elseif (it.eq.1) then         ! b.length

          gdeyrg(iv)= -(ex*xfiv+ey*yfiv+ez*zfiv)

        endif

        if (tesgrd) call gdtreg(nml,iv)

      enddo  ! ... variables in desc. order

      return
      end
! **************************
      subroutine gdtgbl(nml)
!
! CALLS: bldmol,enyreg
!
! -------------------------- gradtest for 'gbpr'

      include 'INCL.H'

      parameter (del=1.d-7)


      ii=(nml-1)*6

      do i = 1,6

! ----------------------------- modify
        pro = gbpr(i,nml)
        gbpr(i,nml) = pro+del
        call bldmol(nml)

        gdn = ( enyreg(nml) - eyrg ) / del

        write (*,*) ' Gb. var #',(ii+i),': ',gdeygb(ii+i),gdn,
     &                                   abs(gdn-gdeygb(ii+i))
! ----------------------------- restore
        gbpr(i,nml) = pro
        call bldmol(nml)

      enddo  ! pars.

      return
      end
! *****************************
      subroutine gdtreg(nml,iv)

! .................................................................
! PURPOSE: calculate partial derivative of reg. energy for molecule
!          'nml' vs. variable 'iv' NUMERICALLY and compare with
!          its value obtained analytically
!
! CALLS:  setvar, enyreg
! .................................................................

      include 'INCL.H'

      parameter (del=1.d-6)

      dimension vlvrx(mxvr)

! ____________________________ get & save values of variables
      do i=1,ivrml1(ntlml)+nvrml(ntlml)-1
        it=ityvr(i)  ! type
        if (it.eq.3) then      ! torsion
          vlvrx(i)=toat(iatvr(i))
        elseif (it.eq.2) then  ! b.angle
          vlvrx(i)=baat(iatvr(i))
        elseif (it.eq.1) then  ! b.length
          vlvrx(i)=blat(iatvr(i))
        endif
      enddo

      ovr=vlvrx(iv)
      vlvrx(iv)=ovr+del       ! change variable 'iv' by 'del'
      call setvar(nml,vlvrx)

      eynw=enyreg(nml)

      gdn=(eynw-eyrg)/del    ! numerical derivative
      gda=gdeyrg(iv)         ! analytical der.

      write (*,'(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
     &       abs(gda-gdn),')'

! _________________________ restore vars
      vlvrx(iv)=ovr
      call setvar(nml,vlvrx)

      return
      end