C------------------------------------------------------------------------------
C Set of diagnostic routines for checking cumulant terms by numerical
C integration. They are not required unless new correlation terms need 
C to be checked.
C------------------------------------------------------------------------------
      subroutine checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,
     & eel_loc_ij)
C Calculate third-order correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.LOCAL'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VECTORS'
      include 'COMMON.FFIELD'
!      real*8 mu(2,maxres),muder(2,maxres),
!     &    muij(4),mu1(2,maxres),mu2(2,maxres),auxvec(2)
       real*8 muij(4),auxvec(2)
      iti=itortyp(itype(i))
      itj=itortyp(itype(j))
      eel_loc_1=a22*b1(1,iti)*b1(1,itj)+a23*b1(1,iti)*b1(2,itj)+
     & a32*b1(2,iti)*b1(1,itj)+a33*b1(2,iti)*b1(2,itj)
      eel_loc_2=a22*b1(1,iti)*Ub2(1,j)+a23*b1(1,iti)*Ub2(2,j)+
     & a32*b1(2,iti)*Ub2(1,j)+a33*b1(2,iti)*Ub2(2,j)
      eel_loc_3=a22*Ub2(1,i)*b1(1,itj)+a23*Ub2(1,i)*b1(2,itj)+
     & a32*Ub2(2,i)*b1(1,itj)+a33*Ub2(2,i)*b1(2,itj)
      eel_loc_4=a22*Ub2(1,i)*Ub2(1,j)+a23*Ub2(1,i)*Ub2(2,j)+
     & a32*Ub2(2,i)*Ub2(1,j)+a33*Ub2(2,i)*Ub2(2,j)
      if (i.gt.iatel_s) then
        iti1=itortyp(itype(i))
      else
        iti1=4
      endif
      iti2=itortyp(itype(i+1))
      itj1=itortyp(itype(j)) 
      if (j.lt.iatel_e+2) then
        itj2=itortyp(itype(j+1))
      else
        itj2=4
      endif
      if (j.lt.nres-1) then
      call integral3(phi(i+2),phi(j+2),iti1,iti2,itj1,itj2,
     & acipa,.false.,eel_1,eel_2,eel_3,eel_4)
      else
      call integral3(phi(i+2),phi(j+2),iti1,iti2,itj1,itj2,
     & acipa,.true.,eel_1,eel_2,eel_3,eel_4)
      endif
cd      write (iout,*) 'eel_1',eel_loc_1,' eel_1_num',4*eel_1
cd      write (iout,*) 'eel_2',eel_loc_2,' eel_2_num',4*eel_2
cd      write (iout,*) 'eel_3',eel_loc_3,' eel_3_num',4*eel_3
cd      write (iout,*) 'eel_4',eel_loc_4,' eel_4_num',4*eel_4
      write (iout,*) 'eel',eel_loc_ij,' eel_num',
     &4*(eel_1+eel_2+eel_3+eel_4)
       return
       end
c----------------------------------------------------------------------
      subroutine checkint4(i,j,k,l,jj,kk,eel4_num)
c Calculate fourth-order correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision gx(3),gx1(3)
      double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
     & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
     & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
      itk = itortyp(itype(k))
C Check integrals
C Copy dipole matrices to temporary arrays
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
          aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
        enddo
      enddo
C Apply inverse transformation
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (j.lt.nres-1) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      if (k.lt.nres-1) then
        do iii=1,2
          aa2(1,iii)=-aa2(1,iii)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      if (l.lt.nres-1) then
        do iii=1,2
          aa2(iii,1)=-aa2(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      if (l.eq.j+1) then
        itl = itortyp(itype(l))
c Compute numerical integrals
        print *,phi(k+2),phi(l+2),itk,itl
        if (l.lt.nres-1) then
cd          write(2,*)'1 ',itk,itl,a_chuj(:,:,jj,i),a_chuj(:,:,kk,k) 
c          call integral(0.0d0,0.0d0,itk,itl,aa1(1,1),
c     &     aa2(1,1),1.0d0,1.0d0,-1.0d0,-1.0d0,.false.,eel4_num)
          call integral(0.0d0,phi(k+2)-pi,0.0d0,phi(l+2)-pi,itk,itl,
     &     aa1(1,1),aa2(1,1),
     &     1.0d0,-1.0d0,1.0d0,-1.0d0,.false.,eel4_num)
        else
cd          write(2,*)'2 ',itk,itl,a_chuj(:,:,jj,i),a_chuj(:,:,kk,k) 
c          call integral(0.0d0,0.0d0,itk,itl,aa1(1,1),
c     &     aa2(1,1),1.0d0,1.0d0,1.0d0,1.0d0,.false.,eel4_num)
          call integral(0.0d0,phi(k+2)-pi,0.0d0,0.0d0,itk,itl,
     &     aa1(1,1),aa2(1,1),
     &     1.0d0,-1.0d0,1.0d0,-1.0d0,.false.,eel4_num)
        endif 
      else
        itl = itortyp(itype(j))
        if (j.lt.nres-1) then
          call integral(0.0d0,phi(k+2)-pi,phi(j+2)-pi,0.0d0,itk,itl,
     &     aa1(1,1),aa2(1,1),1.0d0,-1.0d0,-1.0d0,1.0d0,.true.,eel4_num)
        else
          call integral(0.0d0,phi(k+2)-pi,0.0d0,0.0d0,itk,itl,aa1(1,1),
     &     aa2(1,1),1.0d0,-1.0d0,-1.0d0,1.0d0,.true.,eel4_num)
        endif
      endif
c end check
      return
      end
c-----------------------------------------------------------------------------
      subroutine checkint5(i,j,k,l,jj,kk,eel5_1_num,eel5_2_num,
     &   eel5_3_num,eel5_4_num)
c Calculate fifth-order correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision gx(3),gx1(3)
      double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
     & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
     & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
      iti = itortyp(itype(i))
      itk = itortyp(itype(k))
      itk1= itortyp(itype(k+1))
C Check integrals
C Copy dipole matrices to temporary arrays
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
          aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
        enddo
      enddo
C Apply inverse transformation
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (j.lt.nres-1) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      if (k.lt.nres-1) then
        do iii=1,2
          aa2(1,iii)=-aa2(1,iii)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      if (l.lt.nres-1) then
        do iii=1,2
          aa2(iii,1)=-aa2(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      eel5_1_num=0.0d0
      eel5_2_num=0.0d0
      eel5_3_num=0.0d0
      eel5_4_num=0.0d0
      if (l.eq.j+1) then
        itj = itortyp(itype(j))
        itl = itortyp(itype(l))
        itl1= itortyp(itype(l+1))
c Compute numerical integrals
        if (l.lt.nres-1) then
          if (i.gt.1) then 
            call integral5(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &      1,1,1,1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          else
            call integral5(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &     -1,1,1,1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          endif
        else
          if (i.gt.1) then 
            call integral5(phi(i+2),phi(k+2),phi(j+2),pi,
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &     1,1,1,-1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          else
            call integral5(phi(i+2),phi(k+2),phi(j+2),pi,
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &    -1,1,1,-1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          endif
        endif 
      else
        itj = itortyp(itype(j))
        itl = itortyp(itype(l))
        itj1= itortyp(itype(j+1))
        if (j.lt.nres-1) then
          if (i.gt.1) then 
            call integral5(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &      1,1,1,1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          else
            call integral5(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &     -1,1,1,1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          endif
        else
          if (i.gt.1) then 
            call integral5(phi(i+2),phi(k+2),phi(l+2),pi,
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &     1,1,1,-1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          else
            call integral5(phi(i+2),phi(k+2),phi(l+2),pi,
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &    -1,1,1,-1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
          endif
        endif 
      endif
c end check
      return
      end
c-----------------------------------------------------------------------------
      subroutine checkint6(i,j,k,l,jj,kk,eel6_1_num,eel6_2_num,
     &   eel6_3_num,eel6_4_num,eel6_5_num,eel6_6_num)
c Calculate sixth-order correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision gx(3),gx1(3)
      double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
     & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
     & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
      iti = itortyp(itype(i))
      itk = itortyp(itype(k))
      itk1= itortyp(itype(k+1))
C Check integrals
C Copy dipole matrices to temporary arrays
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
          aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
        enddo
      enddo
C Apply inverse transformation
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (j.lt.nres-1) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      if (k.lt.nres-1) then
        do iii=1,2
          aa2(1,iii)=-aa2(1,iii)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      if (l.lt.nres-1) then
        do iii=1,2
          aa2(iii,1)=-aa2(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      eel6_1_num=0.0d0
      eel6_2_num=0.0d0
      eel6_3_num=0.0d0
      eel6_4_num=0.0d0
      eel6_5_num=0.0d0
      eel6_6_num=0.0d0
      if (l.eq.j+1) then
        itj = itortyp(itype(j))
        itl = itortyp(itype(l))
        itl1= itortyp(itype(l+1))
c Compute numerical integrals
        if (l.lt.nres-1) then
          if (i.gt.1) then 
            call integral6(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &      1,1,1,1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &      eel6_5_num,eel6_6_num)
          else
            call integral6(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &     -1,1,1,1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &      eel6_5_num,eel6_6_num)
          endif
        else
          if (i.gt.1) then 
            call integral6(phi(i+2),phi(k+2),phi(j+2),pi,
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &     1,1,1,-1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &     eel6_5_num,eel6_6_num)
          else
            call integral6(phi(i+2),phi(k+2),phi(j+2),pi,
     &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
     &    -1,1,1,-1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &     eel6_5_num,eel6_6_num)
          endif
        endif 
      else
        itj = itortyp(itype(j))
        itl = itortyp(itype(l))
        itj1= itortyp(itype(j+1))
        if (j.lt.nres-1) then
          if (i.gt.1) then 
            call integral6(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &      1,1,1,1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &      eel6_5_num,eel6_6_num)
          else
            call integral6(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &     -1,1,1,1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &      eel6_5_num,eel6_6_num)
          endif
        else
          if (i.gt.1) then 
            call integral6(phi(i+2),phi(k+2),phi(l+2),pi,
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &     1,1,1,-1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &     eel6_5_num,eel6_6_num)
          else
            call integral6(phi(i+2),phi(k+2),phi(l+2),pi,
     &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
     &    -1,1,1,-1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
     &     eel6_5_num,eel6_6_num)
          endif
        endif 
      endif
c end check
      return
      end
c-----------------------------------------------------------------------------
      subroutine checkint_turn6(i,jj,kk,eel_turn6_num)
c Calculate sixth-order turn correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision gx(3),gx1(3)
      double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
     & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
     & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
      k = i+1
      l = i+3
      j = i+4
      iti = itortyp(itype(i))
      itk = itortyp(itype(k))
      itk1= itortyp(itype(k+1))
C Check integrals
C Copy dipole matrices to temporary arrays
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
          aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
        enddo
      enddo
C Apply inverse transformation
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (j.lt.nres-1) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      if (k.lt.nres-1) then
        do iii=1,2
          aa2(1,iii)=-aa2(1,iii)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      if (l.lt.nres-1) then
        do iii=1,2
          aa2(iii,1)=-aa2(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa2(iii,jjj)=-aa2(iii,jjj)
          enddo
        enddo
      endif
      eel_turn6_num=0.0d0
      itj = itortyp(itype(j))
      itl = itortyp(itype(l))
      itj1= itortyp(itype(j+1))
      call integral_turn6(phi(i+2),phi(i+3),phi(i+4),phi(i+5),
     &  iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),eel_turn6_num)
      write (2,*) 'eel_turn6_num',eel_turn6_num
c end check
      return
      end
c-----------------------------------------------------------------------------
      subroutine checkint_turn3(i,a_temp,eel_turn3_num)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
c Calculate third-order turn correlation terms by numerical integration.
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision a_temp(2,2),aa1(2,2)
      iti1 = itortyp(itype(i+1))
      iti2 = itortyp(itype(i+2))
C Check integrals
C Apply inverse transformation
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_temp(iii,jjj)
        enddo
      enddo
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (i.lt.nres-3) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      eel_turn3_num=0.0d0
c Compute numerical integrals
      if (i.lt.nres-3) then
        call integral3a(phi(i+3),phi(i+4),iti1,iti2,
     &   aa1(1,1), 1,eel_turn3_num)
      else
        call integral3a(phi(i+3),phi(i+4),iti1,iti2,
     &   aa1(1,1),-1,eel_turn3_num)
      endif
c end check
      return
      end
c-----------------------------------------------------------------------------
      subroutine checkint_turn4(i,a_temp,eel_turn4_num)
c Calculate fourth-order turn correlation terms by numerical integration.
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      double precision a_temp(2,2),aa1(2,2)
      iti1 = itortyp(itype(i+1))
      iti2 = itortyp(itype(i+2))
      iti3 = itortyp(itype(i+3))
C Check integrals
C Apply inverse transformation
      do iii=1,2
        do jjj=1,2
          aa1(iii,jjj)=a_temp(iii,jjj)
        enddo
      enddo
      do iii=1,2
        aa1(1,iii)=-aa1(1,iii)
      enddo
      if (i.lt.nres-4) then
        do iii=1,2
          aa1(iii,1)=-aa1(iii,1)
        enddo
      else
        do iii=1,2
          do jjj=1,2
            aa1(iii,jjj)=-aa1(iii,jjj)
          enddo
        enddo
      endif
      eel_turn4_num=0.0d0
c Compute numerical integrals
      if (i.lt.nres-4) then
        call integral4a(phi(i+3),phi(i+4),phi(i+5),
     &   iti1,iti2,iti3,aa1(1,1),1,eel_turn4_num)
      else
        call integral4a(phi(i+3),phi(i+4),phi(i+5),
     &   iti1,iti2,iti3,aa1(1,1),-1,eel_turn4_num)
      endif
c end check
      return
      end