source/unres/src_MD-M/ecorr_num.f

   1 C------------------------------------------------------------------------------
   2 C Set of diagnostic routines for checking cumulant terms by numerical
   3 C integration. They are not required unless new correlation terms need
   4 C to be checked.
   5 C------------------------------------------------------------------------------
   6       subroutine checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,
   7      & eel_loc_ij)
   8 C Calculate third-order correlation terms by numerical integration.
   9       implicit real*8 (a-h,o-z)
  10       include 'DIMENSIONS'
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.GEO'
  13       include 'COMMON.VAR'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.DERIV'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.CONTACTS'
  19       include 'COMMON.TORSION'
  20       include 'COMMON.VECTORS'
  21       include 'COMMON.FFIELD'
  22 !      real*8 mu(2,maxres),muder(2,maxres),
  23 !     &    muij(4),mu1(2,maxres),mu2(2,maxres),auxvec(2)
  24        real*8 muij(4),auxvec(2)
  25       iti=itortyp(itype(i))
  26       itj=itortyp(itype(j))
  27       eel_loc_1=a22*b1(1,iti)*b1(1,itj)+a23*b1(1,iti)*b1(2,itj)+
  28      & a32*b1(2,iti)*b1(1,itj)+a33*b1(2,iti)*b1(2,itj)
  29       eel_loc_2=a22*b1(1,iti)*Ub2(1,j)+a23*b1(1,iti)*Ub2(2,j)+
  30      & a32*b1(2,iti)*Ub2(1,j)+a33*b1(2,iti)*Ub2(2,j)
  31       eel_loc_3=a22*Ub2(1,i)*b1(1,itj)+a23*Ub2(1,i)*b1(2,itj)+
  32      & a32*Ub2(2,i)*b1(1,itj)+a33*Ub2(2,i)*b1(2,itj)
  33       eel_loc_4=a22*Ub2(1,i)*Ub2(1,j)+a23*Ub2(1,i)*Ub2(2,j)+
  34      & a32*Ub2(2,i)*Ub2(1,j)+a33*Ub2(2,i)*Ub2(2,j)
  35       if (i.gt.iatel_s) then
  36         iti1=itortyp(itype(i))
  37       else
  38         iti1=4
  39       endif
  40       iti2=itortyp(itype(i+1))
  41       itj1=itortyp(itype(j))
  42       if (j.lt.iatel_e+2) then
  43         itj2=itortyp(itype(j+1))
  44       else
  45         itj2=4
  46       endif
  47       if (j.lt.nres-1) then
  48       call integral3(phi(i+2),phi(j+2),iti1,iti2,itj1,itj2,
  49      & acipa,.false.,eel_1,eel_2,eel_3,eel_4)
  50       else
  51       call integral3(phi(i+2),phi(j+2),iti1,iti2,itj1,itj2,
  52      & acipa,.true.,eel_1,eel_2,eel_3,eel_4)
  53       endif
  54 cd      write (iout,*) 'eel_1',eel_loc_1,' eel_1_num',4*eel_1
  55 cd      write (iout,*) 'eel_2',eel_loc_2,' eel_2_num',4*eel_2
  56 cd      write (iout,*) 'eel_3',eel_loc_3,' eel_3_num',4*eel_3
  57 cd      write (iout,*) 'eel_4',eel_loc_4,' eel_4_num',4*eel_4
  58       write (iout,*) 'eel',eel_loc_ij,' eel_num',
  59      &4*(eel_1+eel_2+eel_3+eel_4)
  60        return
  61        end
  62 c----------------------------------------------------------------------
  63       subroutine checkint4(i,j,k,l,jj,kk,eel4_num)
  64 c Calculate fourth-order correlation terms by numerical integration.
  65       implicit real*8 (a-h,o-z)
  66       include 'DIMENSIONS'
  67       include 'COMMON.IOUNITS'
  68       include 'COMMON.CHAIN'
  69       include 'COMMON.DERIV'
  70       include 'COMMON.INTERACT'
  71       include 'COMMON.CONTACTS'
  72       include 'COMMON.TORSION'
  73       include 'COMMON.VAR'
  74       include 'COMMON.GEO'
  75       double precision gx(3),gx1(3)
  76       double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
  77      & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
  78      & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
  79       itk = itortyp(itype(k))
  80 C Check integrals
  81 C Copy dipole matrices to temporary arrays
  82       do iii=1,2
  83         do jjj=1,2
  84           aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
  85           aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
  86         enddo
  87       enddo
  88 C Apply inverse transformation
  89       do iii=1,2
  90         aa1(1,iii)=-aa1(1,iii)
  91       enddo
  92       if (j.lt.nres-1) then
  93         do iii=1,2
  94           aa1(iii,1)=-aa1(iii,1)
  95         enddo
  96       else
  97         do iii=1,2
  98           do jjj=1,2
  99             aa1(iii,jjj)=-aa1(iii,jjj)
 100           enddo
 101         enddo
 102       endif
 103       if (k.lt.nres-1) then
 104         do iii=1,2
 105           aa2(1,iii)=-aa2(1,iii)
 106         enddo
 107       else
 108         do iii=1,2
 109           do jjj=1,2
 110             aa2(iii,jjj)=-aa2(iii,jjj)
 111           enddo
 112         enddo
 113       endif
 114       if (l.lt.nres-1) then
 115         do iii=1,2
 116           aa2(iii,1)=-aa2(iii,1)
 117         enddo
 118       else
 119         do iii=1,2
 120           do jjj=1,2
 121             aa2(iii,jjj)=-aa2(iii,jjj)
 122           enddo
 123         enddo
 124       endif
 125       if (l.eq.j+1) then
 126         itl = itortyp(itype(l))
 127 c Compute numerical integrals
 128         print *,phi(k+2),phi(l+2),itk,itl
 129         if (l.lt.nres-1) then
 130 cd          write(2,*)'1 ',itk,itl,a_chuj(:,:,jj,i),a_chuj(:,:,kk,k)
 131 c          call integral(0.0d0,0.0d0,itk,itl,aa1(1,1),
 132 c     &     aa2(1,1),1.0d0,1.0d0,-1.0d0,-1.0d0,.false.,eel4_num)
 133           call integral(0.0d0,phi(k+2)-pi,0.0d0,phi(l+2)-pi,itk,itl,
 134      &     aa1(1,1),aa2(1,1),
 135      &     1.0d0,-1.0d0,1.0d0,-1.0d0,.false.,eel4_num)
 136         else
 137 cd          write(2,*)'2 ',itk,itl,a_chuj(:,:,jj,i),a_chuj(:,:,kk,k)
 138 c          call integral(0.0d0,0.0d0,itk,itl,aa1(1,1),
 139 c     &     aa2(1,1),1.0d0,1.0d0,1.0d0,1.0d0,.false.,eel4_num)
 140           call integral(0.0d0,phi(k+2)-pi,0.0d0,0.0d0,itk,itl,
 141      &     aa1(1,1),aa2(1,1),
 142      &     1.0d0,-1.0d0,1.0d0,-1.0d0,.false.,eel4_num)
 143         endif
 144       else
 145         itl = itortyp(itype(j))
 146         if (j.lt.nres-1) then
 147           call integral(0.0d0,phi(k+2)-pi,phi(j+2)-pi,0.0d0,itk,itl,
 148      &     aa1(1,1),aa2(1,1),1.0d0,-1.0d0,-1.0d0,1.0d0,.true.,eel4_num)
 149         else
 150           call integral(0.0d0,phi(k+2)-pi,0.0d0,0.0d0,itk,itl,aa1(1,1),
 151      &     aa2(1,1),1.0d0,-1.0d0,-1.0d0,1.0d0,.true.,eel4_num)
 152         endif
 153       endif
 154 c end check
 155       return
 156       end
 157 c-----------------------------------------------------------------------------
 158       subroutine checkint5(i,j,k,l,jj,kk,eel5_1_num,eel5_2_num,
 159      &   eel5_3_num,eel5_4_num)
 160 c Calculate fifth-order correlation terms by numerical integration.
 161       implicit real*8 (a-h,o-z)
 162       include 'DIMENSIONS'
 163       include 'COMMON.IOUNITS'
 164       include 'COMMON.CHAIN'
 165       include 'COMMON.DERIV'
 166       include 'COMMON.INTERACT'
 167       include 'COMMON.CONTACTS'
 168       include 'COMMON.TORSION'
 169       include 'COMMON.VAR'
 170       include 'COMMON.GEO'
 171       double precision gx(3),gx1(3)
 172       double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
 173      & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
 174      & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
 175       iti = itortyp(itype(i))
 176       itk = itortyp(itype(k))
 177       itk1= itortyp(itype(k+1))
 178 C Check integrals
 179 C Copy dipole matrices to temporary arrays
 180       do iii=1,2
 181         do jjj=1,2
 182           aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
 183           aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
 184         enddo
 185       enddo
 186 C Apply inverse transformation
 187       do iii=1,2
 188         aa1(1,iii)=-aa1(1,iii)
 189       enddo
 190       if (j.lt.nres-1) then
 191         do iii=1,2
 192           aa1(iii,1)=-aa1(iii,1)
 193         enddo
 194       else
 195         do iii=1,2
 196           do jjj=1,2
 197             aa1(iii,jjj)=-aa1(iii,jjj)
 198           enddo
 199         enddo
 200       endif
 201       if (k.lt.nres-1) then
 202         do iii=1,2
 203           aa2(1,iii)=-aa2(1,iii)
 204         enddo
 205       else
 206         do iii=1,2
 207           do jjj=1,2
 208             aa2(iii,jjj)=-aa2(iii,jjj)
 209           enddo
 210         enddo
 211       endif
 212       if (l.lt.nres-1) then
 213         do iii=1,2
 214           aa2(iii,1)=-aa2(iii,1)
 215         enddo
 216       else
 217         do iii=1,2
 218           do jjj=1,2
 219             aa2(iii,jjj)=-aa2(iii,jjj)
 220           enddo
 221         enddo
 222       endif
 223       eel5_1_num=0.0d0
 224       eel5_2_num=0.0d0
 225       eel5_3_num=0.0d0
 226       eel5_4_num=0.0d0
 227       if (l.eq.j+1) then
 228         itj = itortyp(itype(j))
 229         itl = itortyp(itype(l))
 230         itl1= itortyp(itype(l+1))
 231 c Compute numerical integrals
 232         if (l.lt.nres-1) then
 233           if (i.gt.1) then
 234             call integral5(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
 235      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 236      &      1,1,1,1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 237           else
 238             call integral5(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
 239      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 240      &     -1,1,1,1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 241           endif
 242         else
 243           if (i.gt.1) then
 244             call integral5(phi(i+2),phi(k+2),phi(j+2),pi,
 245      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 246      &     1,1,1,-1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 247           else
 248             call integral5(phi(i+2),phi(k+2),phi(j+2),pi,
 249      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 250      &    -1,1,1,-1,.false.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 251           endif
 252         endif
 253       else
 254         itj = itortyp(itype(j))
 255         itl = itortyp(itype(l))
 256         itj1= itortyp(itype(j+1))
 257         if (j.lt.nres-1) then
 258           if (i.gt.1) then
 259             call integral5(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
 260      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 261      &      1,1,1,1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 262           else
 263             call integral5(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
 264      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 265      &     -1,1,1,1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 266           endif
 267         else
 268           if (i.gt.1) then
 269             call integral5(phi(i+2),phi(k+2),phi(l+2),pi,
 270      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 271      &     1,1,1,-1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 272           else
 273             call integral5(phi(i+2),phi(k+2),phi(l+2),pi,
 274      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 275      &    -1,1,1,-1,.true.,eel5_1_num,eel5_2_num,eel5_3_num,eel5_4_num)
 276           endif
 277         endif
 278       endif
 279 c end check
 280       return
 281       end
 282 c-----------------------------------------------------------------------------
 283       subroutine checkint6(i,j,k,l,jj,kk,eel6_1_num,eel6_2_num,
 284      &   eel6_3_num,eel6_4_num,eel6_5_num,eel6_6_num)
 285 c Calculate sixth-order correlation terms by numerical integration.
 286       implicit real*8 (a-h,o-z)
 287       include 'DIMENSIONS'
 288       include 'COMMON.IOUNITS'
 289       include 'COMMON.CHAIN'
 290       include 'COMMON.DERIV'
 291       include 'COMMON.INTERACT'
 292       include 'COMMON.CONTACTS'
 293       include 'COMMON.TORSION'
 294       include 'COMMON.VAR'
 295       include 'COMMON.GEO'
 296       double precision gx(3),gx1(3)
 297       double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
 298      & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
 299      & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
 300       iti = itortyp(itype(i))
 301       itk = itortyp(itype(k))
 302       itk1= itortyp(itype(k+1))
 303 C Check integrals
 304 C Copy dipole matrices to temporary arrays
 305       do iii=1,2
 306         do jjj=1,2
 307           aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
 308           aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
 309         enddo
 310       enddo
 311 C Apply inverse transformation
 312       do iii=1,2
 313         aa1(1,iii)=-aa1(1,iii)
 314       enddo
 315       if (j.lt.nres-1) then
 316         do iii=1,2
 317           aa1(iii,1)=-aa1(iii,1)
 318         enddo
 319       else
 320         do iii=1,2
 321           do jjj=1,2
 322             aa1(iii,jjj)=-aa1(iii,jjj)
 323           enddo
 324         enddo
 325       endif
 326       if (k.lt.nres-1) then
 327         do iii=1,2
 328           aa2(1,iii)=-aa2(1,iii)
 329         enddo
 330       else
 331         do iii=1,2
 332           do jjj=1,2
 333             aa2(iii,jjj)=-aa2(iii,jjj)
 334           enddo
 335         enddo
 336       endif
 337       if (l.lt.nres-1) then
 338         do iii=1,2
 339           aa2(iii,1)=-aa2(iii,1)
 340         enddo
 341       else
 342         do iii=1,2
 343           do jjj=1,2
 344             aa2(iii,jjj)=-aa2(iii,jjj)
 345           enddo
 346         enddo
 347       endif
 348       eel6_1_num=0.0d0
 349       eel6_2_num=0.0d0
 350       eel6_3_num=0.0d0
 351       eel6_4_num=0.0d0
 352       eel6_5_num=0.0d0
 353       eel6_6_num=0.0d0
 354       if (l.eq.j+1) then
 355         itj = itortyp(itype(j))
 356         itl = itortyp(itype(l))
 357         itl1= itortyp(itype(l+1))
 358 c Compute numerical integrals
 359         if (l.lt.nres-1) then
 360           if (i.gt.1) then
 361             call integral6(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
 362      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 363      &      1,1,1,1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 364      &      eel6_5_num,eel6_6_num)
 365           else
 366             call integral6(phi(i+2),phi(k+2),phi(j+2),phi(l+2),
 367      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 368      &     -1,1,1,1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 369      &      eel6_5_num,eel6_6_num)
 370           endif
 371         else
 372           if (i.gt.1) then
 373             call integral6(phi(i+2),phi(k+2),phi(j+2),pi,
 374      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 375      &     1,1,1,-1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 376      &     eel6_5_num,eel6_6_num)
 377           else
 378             call integral6(phi(i+2),phi(k+2),phi(j+2),pi,
 379      &      iti,itk,itk1,itj,itl,itl1,aa1(1,1),aa2(1,1),
 380      &    -1,1,1,-1,.false.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 381      &     eel6_5_num,eel6_6_num)
 382           endif
 383         endif
 384       else
 385         itj = itortyp(itype(j))
 386         itl = itortyp(itype(l))
 387         itj1= itortyp(itype(j+1))
 388         if (j.lt.nres-1) then
 389           if (i.gt.1) then
 390             call integral6(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
 391      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 392      &      1,1,1,1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 393      &      eel6_5_num,eel6_6_num)
 394           else
 395             call integral6(phi(i+2),phi(k+2),phi(l+2),phi(j+2),
 396      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 397      &     -1,1,1,1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 398      &      eel6_5_num,eel6_6_num)
 399           endif
 400         else
 401           if (i.gt.1) then
 402             call integral6(phi(i+2),phi(k+2),phi(l+2),pi,
 403      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 404      &     1,1,1,-1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 405      &     eel6_5_num,eel6_6_num)
 406           else
 407             call integral6(phi(i+2),phi(k+2),phi(l+2),pi,
 408      &      iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),
 409      &    -1,1,1,-1,.true.,eel6_1_num,eel6_2_num,eel6_3_num,eel6_4_num,
 410      &     eel6_5_num,eel6_6_num)
 411           endif
 412         endif
 413       endif
 414 c end check
 415       return
 416       end
 417 c-----------------------------------------------------------------------------
 418       subroutine checkint_turn6(i,jj,kk,eel_turn6_num)
 419 c Calculate sixth-order turn correlation terms by numerical integration.
 420       implicit real*8 (a-h,o-z)
 421       include 'DIMENSIONS'
 422       include 'COMMON.IOUNITS'
 423       include 'COMMON.CHAIN'
 424       include 'COMMON.DERIV'
 425       include 'COMMON.INTERACT'
 426       include 'COMMON.CONTACTS'
 427       include 'COMMON.TORSION'
 428       include 'COMMON.VAR'
 429       include 'COMMON.GEO'
 430       double precision gx(3),gx1(3)
 431       double precision ee1t(2,2),ee2t(2,2),ee1ta1(2,2),ee2ta2(2,2),
 432      & ee1ta1_der(2,2,3,5),ee2ta2_der(2,2,3,5),aa1(2,2),aa2(2,2),
 433      & aa2t(2,2),uugk(2,2),uugl(2,2),uugj(2,2),pizda(2,2)
 434       k = i+1
 435       l = i+3
 436       j = i+4
 437       iti = itortyp(itype(i))
 438       itk = itortyp(itype(k))
 439       itk1= itortyp(itype(k+1))
 440 C Check integrals
 441 C Copy dipole matrices to temporary arrays
 442       do iii=1,2
 443         do jjj=1,2
 444           aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
 445           aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
 446         enddo
 447       enddo
 448 C Apply inverse transformation
 449       do iii=1,2
 450         aa1(1,iii)=-aa1(1,iii)
 451       enddo
 452       if (j.lt.nres-1) then
 453         do iii=1,2
 454           aa1(iii,1)=-aa1(iii,1)
 455         enddo
 456       else
 457         do iii=1,2
 458           do jjj=1,2
 459             aa1(iii,jjj)=-aa1(iii,jjj)
 460           enddo
 461         enddo
 462       endif
 463       if (k.lt.nres-1) then
 464         do iii=1,2
 465           aa2(1,iii)=-aa2(1,iii)
 466         enddo
 467       else
 468         do iii=1,2
 469           do jjj=1,2
 470             aa2(iii,jjj)=-aa2(iii,jjj)
 471           enddo
 472         enddo
 473       endif
 474       if (l.lt.nres-1) then
 475         do iii=1,2
 476           aa2(iii,1)=-aa2(iii,1)
 477         enddo
 478       else
 479         do iii=1,2
 480           do jjj=1,2
 481             aa2(iii,jjj)=-aa2(iii,jjj)
 482           enddo
 483         enddo
 484       endif
 485       eel_turn6_num=0.0d0
 486       itj = itortyp(itype(j))
 487       itl = itortyp(itype(l))
 488       itj1= itortyp(itype(j+1))
 489       call integral_turn6(phi(i+2),phi(i+3),phi(i+4),phi(i+5),
 490      &  iti,itk,itk1,itl,itj,itj1,aa1(1,1),aa2(1,1),eel_turn6_num)
 491       write (2,*) 'eel_turn6_num',eel_turn6_num
 492 c end check
 493       return
 494       end
 495 c-----------------------------------------------------------------------------
 496       subroutine checkint_turn3(i,a_temp,eel_turn3_num)
 497       implicit real*8 (a-h,o-z)
 498       include 'DIMENSIONS'
 499 c Calculate third-order turn correlation terms by numerical integration.
 500       include 'COMMON.IOUNITS'
 501       include 'COMMON.CHAIN'
 502       include 'COMMON.DERIV'
 503       include 'COMMON.INTERACT'
 504       include 'COMMON.CONTACTS'
 505       include 'COMMON.TORSION'
 506       include 'COMMON.VAR'
 507       include 'COMMON.GEO'
 508       double precision a_temp(2,2),aa1(2,2)
 509       iti1 = itortyp(itype(i+1))
 510       iti2 = itortyp(itype(i+2))
 511 C Check integrals
 512 C Apply inverse transformation
 513       do iii=1,2
 514         do jjj=1,2
 515           aa1(iii,jjj)=a_temp(iii,jjj)
 516         enddo
 517       enddo
 518       do iii=1,2
 519         aa1(1,iii)=-aa1(1,iii)
 520       enddo
 521       if (i.lt.nres-3) then
 522         do iii=1,2
 523           aa1(iii,1)=-aa1(iii,1)
 524         enddo
 525       else
 526         do iii=1,2
 527           do jjj=1,2
 528             aa1(iii,jjj)=-aa1(iii,jjj)
 529           enddo
 530         enddo
 531       endif
 532       eel_turn3_num=0.0d0
 533 c Compute numerical integrals
 534       if (i.lt.nres-3) then
 535         call integral3a(phi(i+3),phi(i+4),iti1,iti2,
 536      &   aa1(1,1), 1,eel_turn3_num)
 537       else
 538         call integral3a(phi(i+3),phi(i+4),iti1,iti2,
 539      &   aa1(1,1),-1,eel_turn3_num)
 540       endif
 541 c end check
 542       return
 543       end
 544 c-----------------------------------------------------------------------------
 545       subroutine checkint_turn4(i,a_temp,eel_turn4_num)
 546 c Calculate fourth-order turn correlation terms by numerical integration.
 547       implicit real*8 (a-h,o-z)
 548       include 'DIMENSIONS'
 549       include 'COMMON.IOUNITS'
 550       include 'COMMON.CHAIN'
 551       include 'COMMON.DERIV'
 552       include 'COMMON.INTERACT'
 553       include 'COMMON.CONTACTS'
 554       include 'COMMON.TORSION'
 555       include 'COMMON.VAR'
 556       include 'COMMON.GEO'
 557       double precision a_temp(2,2),aa1(2,2)
 558       iti1 = itortyp(itype(i+1))
 559       iti2 = itortyp(itype(i+2))
 560       iti3 = itortyp(itype(i+3))
 561 C Check integrals
 562 C Apply inverse transformation
 563       do iii=1,2
 564         do jjj=1,2
 565           aa1(iii,jjj)=a_temp(iii,jjj)
 566         enddo
 567       enddo
 568       do iii=1,2
 569         aa1(1,iii)=-aa1(1,iii)
 570       enddo
 571       if (i.lt.nres-4) then
 572         do iii=1,2
 573           aa1(iii,1)=-aa1(iii,1)
 574         enddo
 575       else
 576         do iii=1,2
 577           do jjj=1,2
 578             aa1(iii,jjj)=-aa1(iii,jjj)
 579           enddo
 580         enddo
 581       endif
 582       eel_turn4_num=0.0d0
 583 c Compute numerical integrals
 584       if (i.lt.nres-4) then
 585         call integral4a(phi(i+3),phi(i+4),phi(i+5),
 586      &   iti1,iti2,iti3,aa1(1,1),1,eel_turn4_num)
 587       else
 588         call integral4a(phi(i+3),phi(i+4),phi(i+5),
 589      &   iti1,iti2,iti3,aa1(1,1),-1,eel_turn4_num)
 590       endif
 591 c end check
 592       return
 593       end