--- /dev/null
+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
projects / unres.git / blobdiff