Merge branch 'devel' into feature-ga
[unres.git] / source / unres / src_MD / src / old_F / ecorr_num.f
diff --git a/source/unres/src_MD/src/old_F/ecorr_num.f b/source/unres/src_MD/src/old_F/ecorr_num.f
deleted file mode 100644 (file)
index 3afecb9..0000000
+++ /dev/null
@@ -1,593 +0,0 @@
-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