! **************************************************************
!
! This file contains the subroutines:  opesol,gdtsol
!
! 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 opesol(nml)

! ......................................................................
! PURPOSE:  derivatives of solvatation energy vs. internal variables for
!           molecule 'nml'
!
!  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:    esolan, gdtsol
! ......................................................................

      include 'INCL.H'
      double precision xfvr, yfvr, zfvr, xfrvr, yfrvr, zfrvr, esolan, dx
      double precision dy, dz, xfat, yfat, zfat, xi, yi, zi, xfrat
      double precision yfrat, zfrat, xb, yb, zb, ex, ey, ez, xfiv, yfiv
      double precision zfiv, xfriv, yfriv, zfriv

      integer ntlvr, nml, ix2, ifivr, ilavr, i, i1s, i1a, io, iv, it, ia
      integer ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, j, i14
      
      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)')
     &           ' opesol> 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'

      do i=ifivr,ilavr  ! variables of 'nml'
        gdeysl(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

      eysl = esolan(nml)

! -------------------------------------------------- f & g for atoms

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

        dx = -gradan(i,1)   ! f = - sigma * dA(r)/dr ( calc. in esolan() )
        dy = -gradan(i,2)
        dz = -gradan(i,3)

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

        xi = xat(i)
        yi = yat(i)
        zi = zat(i)

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

      enddo

      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)

        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.d0
        yfiv=0.d0
        zfiv=0.d0
        xfriv=0.d0
        yfriv=0.d0
        zfriv=0.d0

        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

          gdeysl(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

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

        endif

        if (tesgrd) call gdtsol(nml,iv)

      enddo  ! ... variables in desc. order

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

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

      include 'INCL.H'
      
      double precision del, vlvrx, ovr, eynw, esolan, gda, gdn

      integer i, it, iv, nml


      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=esolan(nml)
      gda=gdeysl(iv)
      gdn=(eynw-eysl)/del    ! numerical derivative

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

! _________________________ restore vars
      vlvrx(iv)=ovr

      call setvar(nml,vlvrx)

      return
      end