double precision function rdif(n,x,g,ider)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      integer n,nmax
      parameter (nmax=max_paropt)
#ifdef MPI
      include 'mpif.h'
      include 'COMMON.MPI'
      integer ierror
      double precision g_(max_paropt),sum_
#endif
      include "COMMON.IOUNITS"
      include "COMMON.PMFDERIV"
      include "COMMON.WEIGHTS"
      include "COMMON.TORSION"
      double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
     &  b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0:3),
     &  eello_turn3,
     &  eello_turn3_b2(3,2),eello_turn3_b1(3,2),eello_turn3_e1(2,2,2),
     &  eello_turn3_e2(2,2,2),eello_turn3_e0m1(3),eello_turn3_e0m2(3),
     &  c2(0:3),em(2,2,2),em2(2,2,2),e0m(3),e0m2(3),jac(nmax)
      double precision v1_kcc_loc(3,3,2),v2_kcc_loc(3,3,2),
     & etoriv1(3,3,2),etoriv2(3,3,2),d1,d2
      double precision eneelloc3(4),eneelloc3b(3,2,2,2,4)
      double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
     & phi,cost1,cost2,sint1,sint2,
     & cosg,sing,cos2g,sin2g,sum,d3,delta,f,fb11,fb21,fb12,fb22,fc11,
     & fc12,fd11,fd12,etori,sint1sint2,cosphi,sinphi,sumvalc,sumvals
      double precision conv /0.01745329251994329576d0/
      double precision r0pp(3)  /5.2739d0,5.4561d0,5.2261d0/
      double precision epspp(3) /1.6027d0,1.4879d0,0.0779d0/
      double precision dCACA(2) /3.784904d0,3.810180d0/
      double precision facp1(4),facpp(4),facturn3,ddist
      logical lder
      double precision x(n),g(n)
      integer ipoint,iipoint,ip,it1,it2,itp1,itp2,i,j,k,l,ii,jj!,irt1,irt2
      integer ider
#ifdef MPI
c nfun.eq.-1 : no broadcast
c n.eq.-1 : stop function evaluation
c nfun.ne.-1 .and. n.ne.-1 : master distribute variables to slaves, all do their share
c of computations
c      if (nfun.ne.-1) then
c      call MPI_Bcast(n,1,MPI_INTEGER,Master,Comm1,ierror)
c      write (2,*) "n",n
c      call flush(iout)
      call MPI_Bcast(ider,1,MPI_INTEGER,Master,Comm1,ierror)
      if (ider.le.0) return
c      call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION,Master,
c     &    Comm1,ierror)
c      endif
c      write (2,*) "x",(x(i),i=1,n)
#endif
      lder=ider.gt.1
c      write (iout,*) "ider",ider," lder",lder
c      call flush(iout)
c      print *,"rdif npoint",npoint
      IF (lder) THEN
c Zero out gradient and hessian
      do i=1,n
        g(i)=0.0d0
      enddo
      ENDIF

      ipoint=0
      ip=0
      b=0.0d0
      sum=0.0d0

c      x(16)=dabs(x(16))
c      wello=x(16)
c-- Following for ala-ala only
c      x(2)=dabs(x(2))
c      wello=x(2)
c      wturn=x(3)
c      x(16)=dabs(x(16))
c      x(17)=dabs(x(17))
c      wello=x(16)
c#ifndef TORUN
c      wello=1.0d0
c#endif
c      wturn=x(17)

      do it1=0,nloctyp-1
      do it2=-nloctyp+1,nloctyp-1
#ifdef DEBUG
      write (iout,*) "it1",it1," it2",it2," e0",e0(it1,it2)
#endif
c following to test with selected types
c      do it1=1,1
c      do it2=-1,-1

      if (it1.lt.nloctyp-1) then
        itp1=1
      else
        itp1=2
      endif
      if (iabs(it2).lt.nloctyp-1) then
        itp2=1
      else
        itp2=2
      endif
      d1=dCACA(itp1)
      d2=dCACA(itp2)
      d3=dCACA(1)
      if (itp1.eq.1) then
        facturn3=wturn_PMF*dsqrt(epspp(1))*4.2d0**3
      else
        facturn3=wturn_PMF*dsqrt(epspp(2))*4.2d0**3
      endif

      ip=ip+1

      delta=180.0d0
c      e0=x(ip)

      facp1(4)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
      if (itp2.eq.1) then
        facp1(2)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
      else
        facp1(2)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
      endif
      if (itp1.eq.1) then
        facp1(3)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
      else
        facp1(3)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
      endif
      if (itp1+itp2.eq.2) then
        facp1(1)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
      else if (itp1+itp2.eq.3) then
        facp1(1)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
      else
        facp1(1)=wello_PMF*dsqrt(epspp(3))*4.2d0**3
      endif
c      print *,"itp1",itp1," itp2",itp2," facp1",facp1

c      irt1=15+20*it1
c      irt2=15+20*iabs(it2)
c-- The following for ALL residues
c      irt1=17+31*it1
c      irt2=17+31*iabs(it2)
c-- The followinig for ala-ala only
c      irt1=3
c      irt2=3
c       print *,"it1",it1," it2",it2
c       print *,irt1,irt2,irt1+31,irt2+31
c----------------------
c      if (FOURIERSPLIT) then
        do i=1,2
          do k=1,3
            b1m(k,i)=bnew1tor(k,i,it2)
            b1m2(k,i)=bnew2tor(k,i,it2)
            b2m1(k,i)=bnew1tor(k,i,it1)
            b2m(k,i)=bnew2tor(k,i,it1)
          enddo
        enddo
        do i=1,2
          do k=1,2
            cm(k,i)=ccnewtor(k,i,it2)
            dm(k,i)=ddnewtor(k,i,it1)
          enddo
        enddo
        do i=1,2
          do j=1,2
            do k=1,2
              em(k,j,i)=eenewtor(k,j,i,it1)
              em2(k,j,i)=eenewtor(k,j,i,it2)
            enddo
          enddo
        enddo
        do k=1,3
          e0m(k)=e0newtor(k,it1)
          e0m2(k)=e0newtor(k,it2)
        enddo
c Invert coefficients for L-chirality
#ifdef DEBUG
      write(iout,*) "it1",it1," it2",it2
      write(iout,*) "b1i (b1m)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
      enddo
      write(iout,*) "b2j (b1m2)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
      enddo
      write(iout,*) "b1j (b2m1)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
      enddo
      write(iout,*) "b2j (b2m)"
      do i=1,3
        write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
      enddo
      write (iout,*) "c"
      do i=1,2
        write (iout,'(2f10.5)') (cm(i,j),j=1,2)
      enddo
      write (iout,*) "d"
      do i=1,2
        write (iout,'(2f10.5)') (dm(i,j),j=1,2)
      enddo
      print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
#endif
      v1_kcc_loc=0.0d0
      v2_kcc_loc=0.0d0
      etoriv1=0.0d0
      etoriv2=0.0d0
c      cm=0.0d0
c      dm=0.0d0
      do k=1,3
        do l=1,3
          v1_kcc_loc(k,l,1)=b2m(l,1)*b1m(k,1)+b2m(l,2)*b1m(k,2)
          v2_kcc_loc(k,l,1)=b2m(l,2)*b1m(k,1)+b2m(l,1)*b1m(k,2)
        enddo
      enddo
      do k=1,2
        do l=1,2
          v1_kcc_loc(k,l,2)=-(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
          v2_kcc_loc(k,l,2)=-(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
        enddo
      enddo
#define CHUJUN
#ifdef CHUJUN
      if (torsion_pmf) then

c      write (iout,*) "it1",it1," it2",it2," etori",e0(it1,it2),
c     &    " mask",mask_e0(it1,it2)
c      write (iout,*)"mask_bnew1 1",it2,
c     & (mask_bnew1tor(j,1,iabs(it2)),j=1,3)
c      write (iout,*)"mask_bnew1 2",it2,
c     & (mask_bnew1tor(j,2,iabs(it2)),j=1,3)
c      write (iout,*)"mask_bnew2 1",it1,(mask_bnew2tor(j,1,it1),j=1,3)
c      write (iout,*)"mask_bnew2 2",it1,(mask_bnew1tor(j,2,it1),j=1,3)
      iipoint=ipoint
      do ii=1,np(it1,it2)
        ipoint=ipoint+1
        i=ipoint
        theta1=zz(1,i)*conv
        theta2=zz(2,i)*conv
        gam=(zz(3,i)+delta-180)*conv
        cost1=dcos(theta1)
        cost2=dcos(theta2)
        sint1=dsin(theta1)
        sint2=dsin(theta2)
        cosg=dcos(gam)
        sing=dsin(gam)
        cos2g=dcos(2*gam)
        sin2g=dsin(2*gam)
#ifdef OLDCALC
        b11=(b1m(1,1)+b1m(2,1)*cost2+b1m(3,1)*cost2*cost2)*sint2
        b12=(b1m(1,2)+b1m(2,2)*cost2+b1m(3,2)*cost2*cost2)*sint2
        b21=(b2m(1,1)+b2m(2,1)*cost1+b2m(3,1)*cost1*cost1)*sint1
        b22=(b2m(1,2)+b2m(2,2)*cost1+b2m(3,2)*cost1*cost1)*sint1
        d11=(dm(1,1)+dm(2,1)*cost1)*sint1*sint1
        d12=(dm(1,2)+dm(2,2)*cost1)*sint1*sint1
        c11=(cm(1,1)+cm(2,1)*cost2)*sint2*sint2
        c12=(cm(1,2)+cm(2,2)*cost2)*sint2*sint2
#ifdef DEBUG
        write (iout,*) "b11",b11," b12",b12," b21",b21," b22",b22
        write (iout,*) "d11",d11," d12",d12," c11",c11," c12",c12
        write (iout,*) "cosg",cosg," cos2g",cos2g," sing",sing,
     &    " sin2g",sin2g
#endif
        f=e0(it1,it2)+(b21*b11+b22*b12)*cosg+(b22*b11+b21*b12)*sing
     &    -(c11*d11-c12*d12)*cos2g-(c12*d11+c11*d12)*sin2g
        fb11=b21*cosg+b22*sing
        fb12=b22*cosg+b21*sing
        fb21=b11*cosg+b12*sing
        fb22=b12*cosg+b11*sing
        fc11=-d11*cos2g-d12*sin2g
        fc12= d12*cos2g-d11*sin2g
        fd11=-c11*cos2g-c12*sin2g
        fd12= c12*cos2g-c11*sin2g
#endif
c check with the formula from unres
        c1(0)=0.0d0
        c2(0)=0.0d0
        c1(1)=1.0d0
        c2(1)=1.0d0
        do j=2,3
          c1(j)=c1(j-1)*cost1
          c2(j)=c2(j-1)*cost2
        enddo
        etori=e0(it1,it2)
        sint1sint2=1.0d0
        do j=1,2
          cosphi=dcos(j*gam)
          sinphi=dsin(j*gam)
          sint1sint2=sint1sint2*sint1*sint2
          sumvalc=0.0d0
          sumvals=0.0d0
          do k=1,3
            do l=1,3
              sumvalc=sumvalc+v1_kcc_loc(l,k,j)*c1(k)*c2(l)
              sumvals=sumvals+v2_kcc_loc(l,k,j)*c1(k)*c2(l)
              etoriv1(l,k,j)=c1(k)*c2(l)*sint1sint2*cosphi
              etoriv2(l,k,j)=c1(k)*c2(l)*sint1sint2*sinphi
            enddo
          enddo
c          if (j.eq.1) print *,"j",j," sumvlac",sumvalc*sint1sint2,
c     &      b21*b11+b22*b12,
c     &      "sumvals",sumvals*sint1sint2,b22*b11+b21*b12
c          if (j.eq.2) print *,"j",j," sumvlac",sumvalc*sint1sint2,
c     &      -(c11*d11-c12*d12),
c     &      "sumvals",sumvals*sint1sint2,-(c12*d11+c11*d12)
          etori=etori+(sumvalc*cosphi+sumvals*sinphi)*sint1sint2
        enddo
c      do k=1,3
c        do l=1,3
c          v1_kcc(k,l,1)=b1m(k,1)*b2m(l,1)+b1m(k,2)*b2m(l,2)
c          v2_kcc(k,l,1)=b1m(k,1)*b2m(l,2)+b1m(k,2)*b2m(l,1)
c        enddo
c      enddo
c      do k=1,2
c        do l=1,2
c          v1_kcc(k,l,2)=-0.25d0*(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
c          v2_kcc(k,l,2)=-0.25d0*(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
c        enddo
c      enddo
c        print *,"i",i," theta",zz(i,1),zz(i,2)," gamma",zz(i,3),
c     &    " f",f," etori",etori
#define CIPUN
#ifdef CIPUN
        yc(1,i)=etori
#ifdef DEBUG
        write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
     &    y(1,i),f,yc(1,i),y(1,i)-yc(1,i)
#endif
        diff(1,i)=y(1,i)-yc(1,i)
        sum=sum+diff(1,i)*diff(1,i)*w(1,i)
c
c Compute derivatives of the energy in torsional parameters
c
        IF (lder) THEN
c
        do k=1,n
          jac(k)=0.0d0
        enddo
c Free term
        if (mask_e0(it1,it2).gt.0) jac(mask_e0(it1,it2))=1.0d0
c Derivatives in b1 of the second residue
c        jj=-1
        do j=1,3
c          jj=jj+2
          jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
          if (jj.gt.0) then
            do k=1,3
              jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
     &        +etoriv2(j,k,1)*b2m(k,2)
            enddo
          endif
          jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
          if (jj.gt.0) then
            do k=1,3
              if (it2.gt.0) then
                jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,2)
     &          +etoriv2(j,k,1)*b2m(k,1)
              else if (it2.lt.0) then
                jac(jj)=jac(jj)-etoriv1(j,k,1)*b2m(k,2)
     &          -etoriv2(j,k,1)*b2m(k,1)
              endif
            enddo
          endif
        enddo
c Derivatives in b2 of the first residue
        do j=1,3
c          jj=jj+2
          jj=iabs(mask_bnew2tor(j,1,it1))
          if (jj.gt.0) then
            do k=1,3
              jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
     &        +etoriv2(k,j,1)*b1m(k,2)
            enddo
          endif
          jj=iabs(mask_bnew2tor(j,2,it1))
          if (jj.gt.0) then
            do k=1,3
              if (it1.gt.0) then
                jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
     &          +etoriv2(k,j,1)*b1m(k,1)
              endif
            enddo
          endif
        enddo
c Derivatives in c
        do j=1,2
c          jj=jj+2
          jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
          if (jj.gt.0) then
            do k=1,2
c              jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,1)
c     &         -0.25d0*etoriv2(j,k,2)*dm(k,2)
              jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,1)
     &         -etoriv2(j,k,2)*dm(k,2)
            enddo
          endif
          jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
          if (jj.gt.0) then
            do k=1,2
              if (it2.gt.0) then
c                jac(jj+1)=jac(jj)+0.25d0*etoriv1(j,k,2)*dm(k,2)
c     &          -0.25d0*etoriv2(j,k,2)*dm(k,1)
                jac(jj)=jac(jj)+etoriv1(j,k,2)*dm(k,2)
     &          -etoriv2(j,k,2)*dm(k,1)
              else if (it2.lt.0) then
c                jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,2)
c     &          +0.25d0*etoriv2(j,k,2)*dm(k,1)
                jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,2)
     &          +etoriv2(j,k,2)*dm(k,1)
              endif
            enddo
          endif
        enddo
c Derivatives in d
        do j=1,2
c          jj=jj+2
          jj=iabs(mask_ddnewtor(j,1,it1))
          if (jj.gt.0) then
            do k=1,2
c              jac(jj)=jac(jj)-0.25d0*etoriv1(k,j,2)*cm(k,1)
c     &         -0.25d0*etoriv2(k,j,2)*cm(k,2) 
              jac(jj)=jac(jj)-etoriv1(k,j,2)*cm(k,1)
     &         -etoriv2(k,j,2)*cm(k,2) 
            enddo
          endif
          jj=iabs(mask_ddnewtor(j,2,it1))
          if (jj.gt.0) then
            do k=1,2
              if (it1.gt.0) then
c                jac(jj)=jac(jj)+0.25d0*etoriv1(k,j,2)*cm(k,2)
c     &          -0.25d0*etoriv2(k,j,2)*cm(k,1)
                jac(jj)=jac(jj)+etoriv1(k,j,2)*cm(k,2)
     &          -etoriv2(k,j,2)*cm(k,1)
              endif
            enddo
          endif
        enddo
c Contributions to the gradient and the hessian
        do k=1,n
          g(k)=g(k)+diff(1,i)*jac(k)*w(1,i)
c          h(k,k)=h(k,k)+jac(k)**2*w(1,i)
c          do l=1,k-1
c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(1,i)
c          enddo
        enddo
#ifdef DEBUG
        print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
     &   " g",g(:n)
#endif
        ENDIF ! lder
#endif
      enddo  ! ii

      endif ! torsion_pmf

      do i=1,2
        do k=1,3
          b1m(k,i)=bnew1(k,i,it2)
          b1m2(k,i)=bnew2(k,i,it2)
          b2m1(k,i)=bnew1(k,i,it1)
          b2m(k,i)=bnew2(k,i,it1)
        enddo
      enddo
      do i=1,2
        do j=1,2
          do k=1,2
            em(k,j,i)=eenew(k,j,i,it1)
            em2(k,j,i)=eenew(k,j,i,it2)
          enddo
        enddo
      enddo
      do k=1,3
        e0m(k)=e0new(k,it1)
        e0m2(k)=e0new(k,it2)
      enddo

      if (turn_pmf) then

      ipoint=iipoint
      do ii=1,np(it1,it2)
        ipoint=ipoint+1
        i=ipoint
        theta1=zz(1,i)*conv
        theta2=zz(2,i)*conv
        gam=(zz(3,i)+delta-180)*conv
#define CYCUN
#ifdef CYCUN
c eturn3 contributions
        call eturn3PMF(theta1,theta2,gam,b2m1,b1m2,em,em2,e0m,e0m2,
     &  facturn3,d1,d2,d3,eello_turn3,eello_turn3_b2,eello_turn3_b1,
     &  eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
        yc(2,i)=eello_turn3
#ifdef DEBUG
        write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
     &    y(2,i),yc(2,i),y(2,i)-yc(2,i)
#endif
        diff(2,i)=y(2,i)-yc(2,i)
        sum=sum+diff(2,i)*diff(2,i)*w(2,i)

        IF (lder) THEN

        jac=0.0d0
        if (mask_turn_PMF.gt.0) jac(mask_turn_PMF)=eello_turn3/wturn_PMF
c Derivatives of in b1 of the first residue
c        jj=-1
        do j=1,3
          jj=mask_bnew1(j,1,it1)
          if (jj.gt.0) then
c            jj=jj+2
            jac(jj)=jac(jj)+eello_turn3_b2(j,1)
          endif
        enddo
        do j=1,3
          jj=mask_bnew1(j,2,it1)
          if (jj.gt.0) then
            if (it1.ne.0) then
              jac(jj)=jac(jj)+eello_turn3_b2(j,2)
            endif
          endif
        enddo
c Derivatives of in b2 of the second residue
        do j=1,3
          jj=mask_bnew2(j,1,iabs(it2))
          if (jj.gt.0) then
            jac(jj)=jac(jj)+eello_turn3_b1(j,1)
          endif
        enddo
        do j=1,3
          jj=mask_bnew2(j,2,iabs(it2))
          if (jj.gt.0) then 
            if (it2.gt.0) then
              jac(jj)=jac(jj)+eello_turn3_b1(j,2)
            else if (it2.lt.0) then
              jac(jj)=jac(jj)-eello_turn3_b1(j,2)
            endif
          endif
        enddo
#define EE
#ifdef EE
c Derivatives in e of the first residue
c        jj=21
        do j=1,2
          jj=mask_eenew(j,1,1,it1)
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,1)
          jj=mask_eenew(j,2,2,it1)
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,2)
          if (it1.ne.0) then
            jj=mask_eenew(j,1,2,it1)
            if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,2)
            jj=mask_eenew(j,2,1,it1)
            if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,1)         
          endif
c          jj=jj+4
        enddo
        jj=mask_e0new(1,it1)
        if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(1)
        if (it1.ne.0) then
          jj=mask_e0new(2,it1)
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(2)         
          jj=mask_e0new(3,it1)
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(3)         
        endif
c Derivatives in e of the second residue
c        jj=21
        do j=1,2
          jj=mask_eenew(j,1,1,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,1)
          jj=mask_eenew(j,2,2,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,2)         
          if (it2.gt.0) then
            jj=mask_eenew(j,1,2,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,2)
            jj=mask_eenew(j,2,1,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,1)
          else if (it2.lt.0) then
            jj=mask_eenew(j,1,2,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,1,2)
            jj=mask_eenew(j,2,1,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,2,1)
          endif
c          jj=jj+4
        enddo
        jj=mask_e0new(1,iabs(it2))
        if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(1)
        if (it2.gt.0) then
          jj=mask_e0new(2,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(2)
          jj=mask_e0new(3,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(3)
        else if (it2.lt.0) then
          jj=mask_e0new(2,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(2)
          jj=mask_e0new(3,iabs(it2))
          if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(3)
        endif
#endif
c Contributions to the gradient and the hessian
        do k=1,n
          g(k)=g(k)+diff(2,i)*jac(k)*w(2,i)
c          h(k,k)=h(k,k)+jac(k)**2*w(2,i)
c          do l=1,k-1
c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(2,i)
c          enddo
        enddo
#ifdef DEBUG
        print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
     &   " g",g(:n)
#endif
c
        ENDIF ! lder
#endif
      enddo

      endif ! turn_pmf
c      nfun=nfun+1
c      rdif=sum
c      return
#endif 
      if (eello_pmf) then

      do ii=1,npel(it1,it2)
        eneelloc3b=0.0d0
        ipoint=ipoint+1
        theta1=zz(1,ipoint)*conv
        theta2=zz(2,ipoint)*conv
        ddist=zz(3,ipoint)**3
        thetap1=zz(4,ipoint)*conv
        thetap2=zz(5,ipoint)*conv
        phi=zz(6,ipoint)*conv
        gam1=zz(7,ipoint)*conv
        gam2=zz(8,ipoint)*conv
#ifdef DEBUG
        print *,"ipoint",ipoint," theta1",theta1," theta2",theta2,
     &   " dd",dd," thetap1",thetap1," thetap2", thetap2," phi",phi,
     &   " gam1",gam1," gam2",gam2
#endif
        do k=1,4
          facpp(k)=facp1(k)/ddist
        enddo
#ifdef DEBUG
        print *,"facp1",facp1," facpp",facpp
#endif
        call eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
     &    phi,b2m1,b2m,b1m,b1m2,facpp,eneelloc3,eneelloc3b,lder)
#ifdef DEBUG
        print *,"eneelloc3",eneelloc3
#endif
        do k=1,4
          yc(k,ipoint)=eneelloc3(k)
          diff(k,ipoint)=y(k,ipoint)-yc(k,ipoint)
          sum=sum+diff(k,ipoint)*diff(k,ipoint)*w(k,ipoint)
        enddo
#ifdef DEBUG
        write (2,'(i7,7f7.1,4(2f10.3,5x))') i,(zz(j,ipoint),j=1,7),
     &    (y(k,ipoint),yc(k,ipoint),k=1,4)
#endif
        IF (lder) THEN
        do l=1,4
c          print *,"ipoint",ipoint," l",l
          do k=1,n
            jac(k)=0.0d0
          enddo
          if (mask_ello_PMF.gt.0)
     &    jac(mask_ello_PMF)=jac(mask_ello_PMF)+eneelloc3(l)/wello_PMF
c#ifndef TORUN
c          jac(16)=0.0d0
c#endif
c          jac(2)=eneelloc3(l)/wello
c Derivatives in b1 of the first residue
c          jj=-1
          do j=1,3
            jj=mask_bnew1(j,1,it1)
c            jj=jj+2
            if (jj.gt.0) jac(jj)=eneelloc3b(j,1,1,1,l)
            jj=mask_bnew1(j,2,it1)
            if (jj.gt.0) then
              if (it1.gt.0) then
                jac(jj)=eneelloc3b(j,2,1,1,l)
              else if (it1.lt.0) then
                jac(jj)=-eneelloc3b(j,2,1,1,l)
              endif
            endif
          enddo
c          print *,"bi1: jac",jac(:n)
c Derivatives in b2 of the first residue
          do j=1,3
c            jj=jj+2
            jj=mask_bnew2(j,1,it1)
            if (jj.gt.0) jac(jj)=eneelloc3b(j,1,2,1,l)
            jj=mask_bnew2(j,2,it1)
            if (jj.gt.0) then
              if (it1.gt.0) then
                jac(jj)=eneelloc3b(j,2,2,1,l)
              else if (it1.lt.0) then
                jac(jj)=-eneelloc3b(j,2,2,1,l)
              endif
            endif
          enddo
c          print *,"bi2: jac",jac(:n)
c Derivatives in b1 of the second residue
c          jj=-1
          do j=1,3
            jj=mask_bnew1(j,1,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,1,2,l)
            jj=mask_bnew1(j,2,iabs(it2))
            if (jj.gt.0) then
              if (it2.gt.0) then
                jac(jj)=jac(jj)+eneelloc3b(j,2,1,2,l)
              else if (it2.lt.0) then
                jac(jj)=jac(jj)-eneelloc3b(j,2,1,2,l)
              endif
            endif
          enddo
c          print *,"bj1: jac",jac(:n)
c Derivatives in b2 of the second residue
          do j=1,3
            jj=mask_bnew2(j,1,iabs(it2))
            if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,2,2,l)
            jj=mask_bnew2(j,2,iabs(it2))
            if (jj.gt.0) then
              if (it2.gt.0) then
                jac(jj)=jac(jj)+eneelloc3b(j,2,2,2,l)
              else if (it2.lt.0) then
                jac(jj)=jac(jj)-eneelloc3b(j,2,2,2,l)
              endif
            endif
          enddo
c          print *,"bj2: jac",jac(:n)
          do k=1,n
            g(k)=g(k)+diff(l,ipoint)*jac(k)*w(l,ipoint)
c            h(k,k)=h(k,k)+jac(k)**2*w(l,ipoint)
c            do ll=1,k-1
c              h(k,ll)=h(k,ll)+jac(k)*jac(ll)*w(l,ipoint)
c            enddo
          enddo
#ifdef DEBUG
          print *,"ipoint",ipoint," jac",jac(:n)," w, diff",w(l,ipoint),
     &     diff(1,ipoint)," g",g(:n)
#endif
        enddo ! l
c Contributions to the gradient and the hessian

        ENDIF
      enddo ! ii

      endif ! eello_PMF

      enddo ! it2
      enddo ! it1

#ifdef MPI
c      if (nfun.ne.-1) then
      sum_=sum
c      write (2,*) "sum",sum
      call flush(2)
      call MPI_Reduce(sum_,sum,1,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
c      endif
c      write (iout,*) "after reducing function lder",lder
c      call flush(iout)
      if (lder) then
        g_=g
        call MPI_Reduce(g_(1),g(1),n,MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,Comm1,ierror)
c      write (iout,*) "after reducing gradient"
c      call flush(iout)
      endif
#endif

c      if (nfun.ne.-1) nfun=nfun+1

      rdif=0.5d0*sum
c      write (iout,*) "Exit rdif"
c      call flush(iout)
c      stop
      return
      end
c----------------------------------------------------------------------
      subroutine eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
     &    phi,bi1,bi2,bj1,bj2,facpp,eneelloc3,eneelloc3b,lder)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "COMMON.IOUNITS"
      double precision theta1,theta2,gam1,gam2,dd,thetap1,thetap2,phi
      double precision bi1(3,2),bi2(3,2),bj1(3,2),bj2(3,2),elpp(2,2),
     &    r0pp(2,2)
      double precision eneelloc3(4),eneelloc3b(3,2,2,2,4),emat1(3,3),
     &    emat2(3,3),emati1(3,3),emati2(3,3),ematj1(3,3),ematj2(3,3)
      double precision mui1(3),mui2(3),muj1(3),muj2(3)
      double precision cost1,cost2,sint1,sint2,costp1,sintp1
      double precision facpp(4)
      logical lder
      common /calcdip/ cost1,cost2,sint1,sint2
      double precision pi /3.14159265358979323844d0/
      integer j,k
      mui1=0.0d0
      mui2=0.0d0
      muj1=0.0d0
      muj2=0.0d0
      cost1=dcos(theta1)
      sint1=dsin(theta1)
      cost2=dcos(theta2)
      sint2=dsin(theta2) 
#ifdef DEBUG
      do j=1,3
      print '(4(a,2f10.5))',"bi1",(bi1(j,k),k=1,2)," bi2",
     &  (bi2(j,k),k=1,2),
     & " bj1",(bj1(j,k),k=1,2)," bj2",(bj2(j,k),k=1,2)
      enddo
      print *,"cost1",cost1," sint1",sint1," cost2",cost2," sint2",sint2
#endif
      do j=1,2
        mui1(j+1)=(bi1(1,j)+bi1(2,j)*cost1+bi1(3,j)*cost1*cost1)*sint1
        muj1(j+1)=(bj1(1,j)+bj1(2,j)*cost2+bj1(3,j)*cost2*cost2)*sint2
        mui2(j+1)=(bi2(1,j)+bi2(2,j)*cost1+bi2(3,j)*cost1*cost1)*sint1
        muj2(j+1)=(bj2(1,j)+bj2(2,j)*cost2+bj2(3,j)*cost2*cost2)*sint2
      enddo
c      mui2(2)=-mui2(2)
c      muj2(2)=-muj2(2)
      mui1(2)=-mui1(2)
      muj1(2)=-muj1(2)
#ifdef DEBUG
      print *,"Before rotation"
      print '(a,3f10.5)',"mui1",mui1
      print '(a,3f10.5)',"muj1",muj1
      print '(a,3f10.5)',"mui2",mui2
      print '(a,3f10.5)',"muj2",muj2
#endif
      sintp1=dsin(thetap1)
      costp1=dcos(thetap1)
c      print *,"costp1",costp1," sintp1",sintp1
      emat1(1,1)=sintp1
      emat1(1,2)=0.0d0
      emat1(1,3)=-costp1
      emat1(2,1)=0.0d0
      emat1(2,2)=1.0d0
      emat1(2,3)=0.0d0
      emat1(3,1)=costp1
      emat1(3,2)=0.0d0
      emat1(3,3)=sintp1
      call coorsys(thetap2,phi,emat2)
#ifdef DEBUG
      print *,"emat1"
      print '(3f10.5)',((emat1(j,k),k=1,3),j=1,3)
      print *,"emat2"
      print '(3f10.5)',((emat2(j,k),k=1,3),j=1,3)
#endif
      call rotmatsub(emat1,gam1,emati1)
      call rotmatsub(emat1,gam1,emati2)
#ifdef DEBUG
      print *,"emati1"
      print '(3f10.5)',((emati1(k,j),k=1,3),j=1,3)
      print *,"emati2"
      print '(3f10.5)',((emati2(k,j),k=1,3),j=1,3)
#endif
      call matvecsub(emati1,mui1)
      call matvecsub(emati2,mui2)
      call rotmatsub(emat2,gam2,ematj1)
      call rotmatsub(emat2,gam2,ematj2)
#ifdef DEBUG
      print *,"ematj1"
      print '(3f10.5)',((ematj1(k,j),k=1,3),j=1,3)
      print *,"ematj2"
      print '(3f10.5)',((ematj2(k,j),k=1,3),j=1,3)
#endif
      call matvecsub(ematj1,muj1)
      call matvecsub(ematj2,muj2)
#ifdef DEBUG
      print *,"After rotation"
      print '(a,3f10.5)',"mui1",mui1
      print '(a,3f10.5)',"muj1",muj1
      print '(a,3f10.5)',"mui2",mui2
      print '(a,3f10.5)',"muj2",muj2
#endif
      call edipole(mui1,muj1,emati1,ematj1,facpp(1),eneelloc3(1),
     &  eneelloc3b(1,1,1,1,1),eneelloc3b(1,1,1,2,1),-1,-1,lder)
      call edipole(mui2,muj1,emati2,ematj1,facpp(2),eneelloc3(2),
     &  eneelloc3b(1,1,2,1,2),eneelloc3b(1,1,1,2,2),1,-1,lder)
      call edipole(mui1,muj2,emati1,ematj2,facpp(3),eneelloc3(3),
     &  eneelloc3b(1,1,1,1,3),eneelloc3b(1,1,2,2,3),-1,1,lder)
      call edipole(mui2,muj2,emati2,ematj2,facpp(4),eneelloc3(4),
     &  eneelloc3b(1,1,2,1,4),eneelloc3b(1,1,2,2,4),1,1,lder)
      return
      end
c-----------------------------------------------------------------------------
      subroutine eturn3PMF(theta1,theta2,gam,b1m,b2m,em1,em2,e0m1,e0m2,
     &  facturn,d1,d2,d3,eello_turn3,eello_turn3_b1,eello_turn3_b2,
     &  eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
      implicit none
      double precision theta1,theta2,gam,cost1,cost2,sint1,sint2,cosg,
     &  sing,sint1sq,sint2sq,aux
      double precision b1m(3,2),b2m(3,2),b1(2),b2(2),em1(2,2,2),
     &  em2(2,2,2),e0m1(3),e0m2(3),e1(2,2),e2(2,2),facturn,emat1(3,3),
     &  emat2(3,3),emat1t(3,3),mu1(2),mu2(2),e2gam(2,2),eello_turn3,
     &  eello_turn3_b2(3,2),eello_turn3_b1(3,2),eello_turn3_e1(2,2,2),
     &  eello_turn3_e2(2,2,2),eello_turn3_e0m1(3),eello_turn3_e0m2(3),
     &  r1(3),r2(3),er(3),a(3,3),atemp(3,3),etemp1(2,2),etemp2(2,2),
     &  d1,d2,d3,ddist,dd3
      double precision pi /3.14159265358979323844d0/
      integer i,j,k,l
c diagnostics
c      d1=3.8
c      d2=3.8
c      d3=3.8
c end diagnostics
      cost1=dcos(theta1)
      cost2=dcos(theta2)
      sint1=dsin(theta1)
      sint1sq=sint1*sint1
      sint2=dsin(theta2)
      sint2sq=sint2*sint2
      cosg=dcos(gam)
      sing=dsin(gam)
      mu1(1)=(b1m(1,1)+b1m(2,1)*cost1+b1m(3,1)*cost1*cost1)*sint1
      mu1(2)=(b1m(1,2)+b1m(2,2)*cost1+b1m(3,2)*cost1*cost1)*sint1
      mu2(1)=(b2m(1,1)+b2m(2,1)*cost2+b2m(3,1)*cost2*cost2)*sint2
      mu2(2)=(b2m(1,2)+b2m(2,2)*cost2+b2m(3,2)*cost2*cost2)*sint2
      mu1(1)=-mu1(1)
      mu2(1)=-mu2(1)
      do k=1,2
        do l=1,2
          e1(k,l)=(em1(1,k,l)+em1(2,k,l)*cost1)*sint1sq
          e2(k,l)=(em2(1,k,l)+em2(2,k,l)*cost2)*sint2sq
        enddo
      enddo
      e1(1,1)=e1(1,1)+e0m1(1)*cost1
      e1(1,2)=e1(1,2)+e0m1(2)+e0m1(3)*cost1
      e1(2,1)=e1(2,1)+e0m1(2)*cost1+e0m1(3)
      e1(2,2)=e1(2,2)-e0m1(1)
      e2(1,1)=e2(1,1)+e0m2(1)*cost2
      e2(1,2)=e2(1,2)+e0m2(2)+e0m2(3)*cost2
      e2(2,1)=e2(2,1)+e0m2(2)*cost2+e0m2(3)
      e2(2,2)=e2(2,2)-e0m2(1)
#ifdef DEBUG
      print *,"mu1",mu1," mu2",mu2," e1",e1," e2",e2
#endif
c Coordinate system of the first  dipole
      emat1(1,1)=sint1
      emat1(1,2)=cost1
      emat1(1,3)=0.0d0
      emat1(2,1)=-cost1
      emat1(2,2)=sint1
      emat1(2,3)=0.0d0
      emat1(3,1)=0.0d0
      emat1(3,2)=0.0d0
      emat1(3,3)=1.0d0
      do i=1,3
        do j=1,3
          emat1t(i,j)=emat1(j,i)
        enddo
      enddo
#ifdef DEBUG
      print *,"theta1",theta1*180/pi," theta2",theta2*180/pi,
     &  " gam",gam*180/pi
#endif
c Calculate the reference system of the second unit
      call coorsys1(theta2,gam,emat2)
c Dipole-dipole interaction matrix
#ifdef DEBUG
      print *,"d1",d1," d2",d2," d3",d3
      print *,"emat1",emat1," emat2",emat2
#endif
      r1(1)=-d1/2
      r1(2)=0.0d0
      r1(3)=0.0d0
      call matvecsub(emat1,r1)
      r2(1)=-d3/2
      r2(2)=0.0d0
      r2(3)=0.0d0
#ifdef DEBUG
      print *," r2",r2
#endif
      call matvecsub(emat2,r2)
#ifdef DEBUG
      print *," r2",r2
#endif
      r2(2)=r2(2)+d2
      ddist=dsqrt((r2(1)-r1(1))**2+(r2(2)-r1(2))**2+(r2(3)-r1(3))**2)
      dd3=facturn/ddist**3
      er(1)=(r2(1)-r1(1))/ddist
      er(2)=(r2(2)-r1(2))/ddist
      er(3)=(r2(3)-r1(3))/ddist
#ifdef DEBUG
      print *,"norm",er(1)**2+er(2)**2+er(3)**2
      print *,"r1",r1," r2",r2
      print *," er",er
      print *," dd",dd," dd3",dd3
#endif
      do i=1,3
        do j=1,3
          a(i,j)=-3*er(i)*er(j)*dd3
c          a(i,j)=-3*er(i)*er(j)
        enddo
      enddo
      do i=1,3
        a(i,i)=dd3+a(i,i)
c        a(i,i)=1.0d0+a(i,i)
      enddo
#ifdef DEBUG
      print *,"a before transformation",a
      print *,"emat1t",emat1t
      print *,"emat2",emat2
#endif
      call matmat(emat1t,a,atemp)
#ifdef DEBUG
      print *,"a after transformation",atemp
#endif
      call matmat(atemp,emat2,a)
#ifdef DEBUG
      print *,"a after transformation",a
      print *,"a and mu1"
      print "(1h[,2f10.5,1h],5x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
     &  mu1(1),mu1(2),a(2,2),a(2,3),mu2(1)
      print "(27x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
     &  a(3,2),a(3,3),mu2(2)
#endif
      eello_turn3=(a(2,2)*mu1(1)*mu2(1)+a(2,3)*mu1(1)*mu2(2)
     &            +a(3,2)*mu1(2)*mu2(1)+a(3,3)*mu1(2)*mu2(2))
#ifdef DEBUG
      print *,"eello_turn3_1",eello_turn3
#endif
      eello_turn3_b1(1,1)=-(a(2,2)*mu2(1)+a(2,3)*mu2(2))*sint1
      eello_turn3_b1(2,1)=eello_turn3_b1(1,1)*cost1
      eello_turn3_b1(3,1)=eello_turn3_b1(2,1)*cost1
      eello_turn3_b1(1,2)=(a(3,2)*mu2(1)+a(3,3)*mu2(2))*sint1
      eello_turn3_b1(2,2)=eello_turn3_b1(1,2)*cost1
      eello_turn3_b1(3,2)=eello_turn3_b1(2,2)*cost1
      eello_turn3_b2(1,1)=-(a(2,2)*mu1(1)+a(3,2)*mu1(2))*sint2
      eello_turn3_b2(2,1)=eello_turn3_b2(1,1)*cost2
      eello_turn3_b2(3,1)=eello_turn3_b2(2,1)*cost2
      eello_turn3_b2(1,2)=(a(2,3)*mu1(1)+a(3,3)*mu1(2))*sint2
      eello_turn3_b2(2,2)=eello_turn3_b2(1,2)*cost2
      eello_turn3_b2(3,2)=eello_turn3_b2(2,2)*cost2
      e2gam(1,1)=e2(1,1)*cosg+e2(2,1)*sing
      e2gam(1,2)=e2(1,2)*cosg+e2(2,2)*sing
      e2gam(2,1)=-e2(1,1)*sing+e2(2,1)*cosg
      e2gam(2,2)=-e2(1,2)*sing+e2(2,2)*cosg
      do i=1,2
        do j=1,2
          e2(i,j)=e2gam(i,j)
        enddo
      enddo
c      e1(1,1)=-e1(1,1)
c      e1(1,2)=-e1(1,2)
      e2(1,1)=-e2(1,1)
      e2(1,2)=-e2(1,2)
      a(2,3)=-a(2,3)
      a(3,2)=-a(3,2)
      aux=0.0d0
      do i=1,2
        do j=1,2
          aux=aux+a(i+1,j+1)*(e1(i,1)*e2(1,j)+e1(i,2)*e2(2,j))
        enddo
      enddo
#ifdef DEBUG
      print *,"eello_turn3_2",aux/2
#endif
      eello_turn3=eello_turn3+aux/2 
c      stop
c Derivatives in the elements of complete e1 and e2 matrices
      do i=1,2
        do j=1,2
          etemp1(i,j)=(a(i+1,2)*e2(j,1)+a(i+1,3)*e2(j,2))/2
          etemp2(i,j)=(a(2,j+1)*e1(1,i)+a(3,j+1)*e1(2,i))/2
        enddo
      enddo
c Derivatives in the components of e1
      do i=1,2
        eello_turn3_e1(1,1,i)=etemp1(1,i)*sint1sq
        eello_turn3_e1(2,1,i)=eello_turn3_e1(1,1,i)*cost1
        eello_turn3_e1(1,2,i)= etemp1(2,i)*sint1sq
        eello_turn3_e1(2,2,i)=eello_turn3_e1(1,2,i)*cost1
      enddo
c Derivatives in the components of e01m
      eello_turn3_e0m1(1)=etemp1(1,1)*cost1-etemp1(2,2)
      eello_turn3_e0m1(2)=etemp1(1,2)+cost1*etemp1(2,1)
      eello_turn3_e0m1(3)=etemp1(1,2)*cost1+etemp1(2,1)
c Derivatives in the elements of e2
      do i=1,2
       eello_turn3_e2(1,1,i)=-(etemp2(1,i)*cosg+etemp2(2,i)*sing)
     &   *sint2sq
       eello_turn3_e2(2,1,i)=eello_turn3_e2(1,1,i)*cost2
       eello_turn3_e2(1,2,i)=(-etemp2(1,i)*sing+etemp2(2,i)*cosg)
     &   *sint2sq
       eello_turn3_e2(2,2,i)=eello_turn3_e2(1,2,i)*cost2
      enddo
c Derivatives in the elements of e02m
      eello_turn3_e0m2(1)=-(etemp2(1,1)*cosg+etemp2(2,1)*sing)*cost2
     &                 +etemp2(1,2)*sing-etemp2(2,2)*cosg
      eello_turn3_e0m2(2)=-etemp2(1,2)*cosg-etemp2(2,2)*sing+
     &  (-etemp2(1,1)*sing+etemp2(2,1)*cosg)*cost2
      eello_turn3_e0m2(3)=(-etemp2(1,2)*cosg-etemp2(2,2)*sing)*cost2
     &  -etemp2(1,1)*sing+etemp2(2,1)*cosg
      return
      end
c-----------------------------------------------------------------------------
      subroutine edipole(mui,muj,emati,ematj,facpp,ene,enebi,enebj,
     &  idiri,idirj,lder)
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "COMMON.IOUNITS"
      double precision mui(3),muj(3),emati(3,3),ematj(3,3),
     &  facpp,ene,enebi(3,2),enebj(3,2),muimuj,cost1,sint1,cost2,sint2
      integer idiri,idirj
      logical lder
      common /calcdip/ cost1,cost2,sint1,sint2
      muimuj=mui(1)*muj(1)+mui(2)*muj(2)+mui(3)*muj(3)
      ene=facpp*(muimuj-3*mui(3)*muj(3))
      if (idiri.gt.0) then
        enebi(1,1)=(emati(1,2)*muj(1)+emati(2,2)*muj(2)
     &   +emati(3,2)*muj(3)-3*muj(3)*emati(3,2))*sint1*facpp
        enebi(2,1)=enebi(1,1)*cost1
        enebi(3,1)=enebi(2,1)*cost1
        enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
     &   +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
        enebi(2,2)=enebi(1,2)*cost1
        enebi(3,2)=enebi(2,2)*cost1
      else
        enebi(1,1)=(-emati(1,2)*muj(1)-emati(2,2)*muj(2)
     &   -emati(3,2)*muj(3)+3*muj(3)*emati(3,2))*sint1*facpp
        enebi(2,1)=enebi(1,1)*cost1
        enebi(3,1)=enebi(2,1)*cost1
        enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
     &   +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
        enebi(2,2)=enebi(1,2)*cost1
        enebi(3,2)=enebi(2,2)*cost1
      endif
      if (idirj.gt.0) then
        enebj(1,1)=(ematj(1,2)*mui(1)+ematj(2,2)*mui(2)
     &   +ematj(3,2)*mui(3)-3*mui(3)*ematj(3,2))*sint2*facpp
        enebj(2,1)=enebj(1,1)*cost2
        enebj(3,1)=enebj(2,1)*cost2
        enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
     &   +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
        enebj(2,2)=enebj(1,2)*cost2
        enebj(3,2)=enebj(2,2)*cost2
      else
        enebj(1,1)=(-ematj(1,2)*mui(1)-ematj(2,2)*mui(2)
     &   -ematj(3,2)*mui(3)+3*mui(3)*ematj(3,2))*sint2*facpp
        enebj(2,1)=enebj(1,1)*cost2
        enebj(3,1)=enebj(2,1)*cost2
        enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
     &   +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
        enebj(2,2)=enebj(1,2)*cost2
        enebj(3,2)=enebj(2,2)*cost2
      endif
      return
      end
c-----------------------------------------------------------------------------
      subroutine matvecsub(mat,vec)
      implicit none
      double precision mat(3,3),vec(3),auxvec(3)
      integer i,j
      do i=1,3
        auxvec(i)=0.0d0
        do j=1,3
          auxvec(i)=auxvec(i)+mat(i,j)*vec(j)
        enddo
      enddo
      vec=auxvec
      return
      end
c-------------------------------------------------------------------------------
      subroutine rotmatsub(mat,gam,matgam)
      implicit none
      double precision mat(3,3),gam,matgam(3,3),cosg,sing
      integer i
      cosg=dcos(gam)
      sing=dsin(gam)
      do i=1,3
        matgam(i,1)=mat(i,1)
        matgam(i,2)= mat(i,2)*cosg+mat(i,3)*sing
        matgam(i,3)=-mat(i,2)*sing+mat(i,3)*cosg
      enddo
      return
      end
c-----------------------------------------------------------------------------
      subroutine coorsys(theta,phi,emat)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "COMMON.IOUNITS"
      double precision theta,phi,emat(3,3)
      double precision cost,sint,cosphi,sinphi
c      print *,"theta",theta," phi",phi
      cost=dcos(theta)
      sint=dsin(theta)
      cosphi=dcos(phi)
      sinphi=dsin(phi)
#ifdef DEBUG
      print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
#endif
      emat=0.0d0
      emat(1,1)=sint*cosphi
      emat(2,1)=sint*sinphi
      emat(3,1)=cost
      emat(1,2)=-sinphi
      emat(2,2)=cosphi
      emat(3,2)=0.0d0
      emat(1,3)=-cost*cosphi
      emat(2,3)=-cost*sinphi
      emat(3,3)=sint
      return
      end
c-----------------------------------------------------------------------------
      subroutine coorsys1(theta,phi,emat)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "COMMON.IOUNITS"
      double precision theta,phi,emat(3,3)
      double precision cost,sint,cosphi,sinphi
c      print *,"theta",theta," phi",phi
      cost=dcos(theta)
      sint=dsin(theta)
      cosphi=dcos(phi)
      sinphi=dsin(phi)
#ifdef DEBUG
      print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
#endif
      emat=0.0d0
      emat(1,1)=sint*cosphi
      emat(2,1)=cost
      emat(3,1)=-sint*sinphi
      emat(1,2)=-cost*cosphi
      emat(2,2)=sint
      emat(3,2)=cost*sinphi
      emat(1,3)=sinphi
      emat(2,3)=0.0d0
      emat(3,3)=cosphi
      return
      end
c-------------------------------------------------------------------------------
      subroutine matmat(a,b,c)
      implicit none
      double precision a(3,3),b(3,3),c(3,3)
      integer i,j,k
      do i=1,3
        do j=1,3
          c(i,j)=0.0d0
          do k=1,3
            c(i,j)=c(i,j)+a(i,k)*b(k,j)
          enddo
        enddo
      enddo
      return
      end