-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