source/maxlik/src_FPy.org/angnorm.f

   1       subroutine add_angpair(ici,icj,nang_pair,iang_pair)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'COMMON.IOUNITS'
   6       include 'COMMON.CHAIN'
   7       integer ici,icj,nang_pair,iang_pair(2,maxres)
   8       integer i,ian1,ian2
   9 c      write (iout,*) "add_angpair: ici",ici," icj",icj,
  10 c     &  " nang_pair",nang_pair
  11       ian1=ici+2
  12       if (ian1.lt.4 .or. ian1.gt.nres) return
  13       ian2=icj+2
  14 c      write (iout,*) "ian1",ian1," ian2",ian2
  15       if (ian2.lt.4 .or. ian2.gt.nres) return
  16       do i=1,nang_pair
  17         if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
  18       enddo
  19       nang_pair=nang_pair+1
  20       iang_pair(1,nang_pair)=ian1
  21       iang_pair(2,nang_pair)=ian2
  22       return
  23       end
  24 c-------------------------------------------------------------------------
  25       subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,
  26      &  ib,iprot,lprn)
  27       implicit none
  28       include 'DIMENSIONS'
  29       include 'DIMENSIONS.ZSCOPT'
  30       include 'DIMENSIONS.COMPAR'
  31       include 'COMMON.IOUNITS'
  32       include 'COMMON.VAR'
  33       include 'COMMON.COMPAR'
  34       include 'COMMON.CHAIN'
  35       include 'COMMON.GEO'
  36       integer jfrag,ishif1,ishif2,ib,iprot
  37       integer nn,npart,nn4,nne,i,nbeg,nend,longest,ll,j
  38       double precision deltang,diffang_max,angn,fract
  39       double precision ff
  40       double precision pinorm,spherang
  41       logical lprn
  42       if (lprn) write (iout,'(80(1h*))')
  43       angn=0.0d0
  44       nn = 0
  45       fract = 1.0d0
  46       npart = npiece(jfrag,1,ib,iprot)
  47       nn4 = nstart_sup(iprot)+3
  48       nne = min0(nend_sup(iprot),nres)
  49       if (lprn) write (iout,*) "nn4",nn4," nne",nne
  50       do i=1,npart
  51         nbeg = ifrag(1,i,jfrag,ib,iprot) + 3 - ishif1
  52         if (nbeg.lt.nn4) nbeg=nn4
  53         nend = ifrag(2,i,jfrag,ib,iprot) + 1 - ishif2
  54         if (nend.gt.nne) nend=nne
  55         if (nend.ge.nbeg) then
  56         nn = nn + nend - nbeg + 1
  57         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
  58      &    " nn",nn," ishift1",ishif1," ishift2",ishif2
  59         if (lprn) write (iout,*) "angles"
  60         longest=0
  61         ll = 0
  62         do j=nbeg,nend
  63 c          deltang = pinorm(phi(j)-phi_ref(j+ishif1,iprot))
  64           deltang=spherang(phi_ref(j+ishif1,iprot),
  65      &      theta_ref(j-1+ishif1,iprot),theta_ref(j+ishif1,iprot),
  66      &      phi(j),theta(j-1),theta(j))
  67           if (dabs(deltang).gt.diffang_max) then
  68             if (ll.gt.longest) longest = ll
  69             ll = 0
  70           else
  71             ll=ll+1
  72           endif
  73           if (ll.gt.longest) longest = ll
  74           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
  75      &     rad2deg*phi_ref(j+ishif1,iprot),rad2deg*deltang
  76           angn=angn+dabs(deltang)
  77         enddo
  78         longest=longest+3
  79         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
  80         if (lprn) write (iout,*)"segment",i," longest fragment within",
  81      &    diffang_max*rad2deg,":",longest," fraction",ff
  82         if (ff.lt.fract) fract = ff
  83         endif
  84       enddo
  85       if (nn.gt.0) then
  86         angn = angn/nn
  87       else
  88         angn = dwapi
  89       endif
  90       if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,
  91      &  " fract",fract
  92       return
  93       end
  94 c-------------------------------------------------------------------------
  95       subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,
  96      &  diffang_max,anorm,fract,ib,iprot)
  97       implicit none
  98       include 'DIMENSIONS'
  99       include 'DIMENSIONS.ZSCOPT'
 100       include 'DIMENSIONS.COMPAR'
 101       include 'COMMON.IOUNITS'
 102       include 'COMMON.VAR'
 103       include 'COMMON.COMPAR'
 104       include 'COMMON.CHAIN'
 105       include 'COMMON.GEO'
 106       integer ib,iprot
 107       integer ncont,icont(2,ncont),longest,jfrag,ishif1,ishif2
 108       integer npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
 109       integer nn4,nne,npart,i,jstart,ic1,jend,ic2,idi,iic,j,
 110      &  nn,ishifc,nbeg,nend,ll
 111       double precision diffang_max,anorm,fract
 112       double precision angn, deltang,ff
 113       double precision pinorm,spherang
 114       logical lprn
 115       if (lprn) write (iout,'(80(1h*))')
 116 c
 117 c Determine the segments for which angles will be compared
 118 c
 119       nn4 = nstart_sup(iprot)+3
 120       nne = min0(nend_sup(iprot),nres)
 121       if (lprn) write (iout,*) "nn4",nn4," nne",nne
 122       npart=npiece(jfrag,1,ib,iprot)
 123       npiece_c=0
 124       do i=1,npart
 125 c        write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag,ib,iprot),
 126 c     &    ifrag(2,i,jfrag,ib,iprot)
 127         if (icont(1,ncont).lt.ifrag(1,i,jfrag,ib,iprot) .or.
 128      &    icont(1,1).gt.ifrag(2,i,jfrag,ib,iprot)) goto 11
 129         jstart=1
 130         do while (jstart.lt.ncont .and.
 131      &   icont(1,jstart).lt.ifrag(1,i,jfrag,ib,iprot))
 132 c          write (iout,*) "jstart",jstart," icont",icont(1,jstart),
 133 c     &     " ifrag",ifrag(1,i,jfrag,ib,iprot)
 134           jstart=jstart+1
 135         enddo
 136 c        write (iout,*) "jstart",jstart," icont",icont(1,jstart),
 137 c     &   " ifrag",ifrag(1,i,jfrag,ib,iprot)
 138         if (icont(1,jstart).lt.ifrag(1,i,jfrag,ib,iprot)) goto 11
 139         npiece_c=npiece_c+1
 140         ic1=icont(1,jstart)
 141         ifrag_c(1,npiece_c)=icont(1,jstart)
 142         jend=ncont
 143         do while(jend.gt.1.and.icont(1,jend).gt.
 144      &    ifrag(2,i,jfrag,ib,iprot))
 145 c          write (iout,*) "jend",jend," icont",icont(1,jend),
 146 c     &     " ifrag",ifrag(2,i,jfrag,ib,iprot)
 147           jend=jend-1
 148         enddo
 149 c        write (iout,*) "jend",jend," icont",icont(1,jend),
 150 c     &   " ifrag",ifrag(2,i,jfrag,ib,iprot)
 151         ic2=icont(1,jend)
 152         ifrag_c(2,npiece_c)=icont(1,jend)+1
 153         ishift_c(npiece_c)=ishif1
 154 c        write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
 155 c     &    " ic1",ic1," ic2",ic2,
 156 c     &    " ifrag",ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot)
 157    11   continue
 158         if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
 159           idi=1
 160         else
 161           idi=-1
 162         endif
 163 c        write (iout,*) "idi",idi
 164         if (idi.eq.1) then
 165           if (icont(2,1).gt.ifrag(2,i,jfrag,ib,iprot) .or.
 166      &      icont(2,ncont).lt.ifrag(1,i,jfrag,ib,iprot)) goto 12
 167           jstart=1
 168           do while (jstart.lt.ncont .and.
 169      &     icont(2,jstart).lt.ifrag(1,i,jfrag,ib,iprot))
 170 c           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 171 c     &     " ifrag",ifrag(1,i,jfrag,ib,iprot)
 172             jstart=jstart+1
 173           enddo
 174 c          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 175 c     &     " ifrag",ifrag(1,i,jfrag,ib,iprot)
 176           if (icont(2,jstart).lt.ifrag(1,i,jfrag,ib,iprot)) goto 12
 177           npiece_c=npiece_c+1
 178           ic1=icont(2,jstart)
 179           ifrag_c(2,npiece_c)=icont(2,jstart)+1
 180           jend=ncont
 181           do while(jend.gt.1.and.icont(2,jend).gt.
 182      &      ifrag(2,i,jfrag,ib,iprot))
 183 c            write (iout,*) "jend",jend," icont",icont(2,jend),
 184 c     &     " ifrag",ifrag(2,i,jfrag,ib,iprot)
 185             jend=jend-1
 186           enddo
 187 c          write (iout,*) "jend",jend," icont",icont(2,jend),
 188 c     &     " ifrag",ifrag(2,i,jfrag,ib,iprot)
 189         else if (idi.eq.-1) then
 190           if (icont(2,ncont).gt.ifrag(2,i,jfrag,ib,iprot) .or.
 191      &        icont(2,1).lt.ifrag(1,i,jfrag,ib,iprot)) goto 12
 192           jstart=ncont
 193           do while (jstart.gt.ncont .and.
 194      &     icont(2,jstart).lt.ifrag(1,i,jfrag,ib,iprot))
 195 c           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 196 c     &     " ifrag",ifrag(1,i,jfrag,ib,iprot)
 197             jstart=jstart-1
 198           enddo
 199 c          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 200 c     &     " ifrag",ifrag(1,i,jfrag,ib,iprot)
 201           if (icont(2,jstart).lt.ifrag(1,i,jfrag,ib,iprot)) goto 12
 202           npiece_c=npiece_c+1
 203           ic1=icont(2,jstart)
 204           ifrag_c(2,npiece_c)=icont(2,jstart)+1
 205           jend=1
 206           do while (jend.lt.ncont .and.
 207      &       icont(2,jend).gt.ifrag(2,i,jfrag,ib,iprot))
 208 c             write (iout,*) "jend",jend," icont",icont(2,jend),
 209 c     &         " ifrag",ifrag(2,i,jfrag,ib,iprot)
 210             jend=jend+1
 211           enddo
 212 c          write (iout,*) "jend",jend," icont",icont(2,jend),
 213 c     &     " ifrag",ifrag(2,i,jfrag,ib,iprot)
 214         endif
 215         ic2=icont(2,jend)
 216         if (ic2.lt.ic1) then
 217           iic = ic1
 218           ic1 = ic2
 219           ic2 = iic
 220         endif
 221 c        write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
 222 c     &    " jstart:",jstart," jend",jend,
 223 c     &    " ifrag",ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot)
 224         ifrag_c(1,npiece_c)=ic1
 225         ifrag_c(2,npiece_c)=ic2+1
 226         ishift_c(npiece_c)=ishif2
 227    12   continue
 228       enddo
 229       if (lprn) then
 230         write (iout,*) "Before merge: npiece_c",npiece_c
 231         do i=1,npiece_c
 232           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
 233         enddo
 234       endif
 235 c
 236 c Merge overlapping segments (e.g., avoid splitting helices)
 237 c
 238       i=1
 239       do while (i .lt. npiece_c)
 240         if (ishift_c(i).eq.ishift_c(i+1) .and.
 241      &     ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
 242            ifrag_c(2,i)=ifrag_c(2,i+1)
 243            do j=i+1,npiece_c
 244              ishift_c(j)=ishift_c(j+1)
 245              ifrag_c(1,j)=ifrag_c(1,j+1)
 246              ifrag_c(2,j)=ifrag_c(2,j+1)
 247            enddo
 248            npiece_c=npiece_c-1
 249         else
 250           i=i+1
 251         endif
 252       enddo
 253       if (lprn) then
 254         write (iout,*) "After merge: npiece_c",npiece_c
 255         do i=1,npiece_c
 256           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
 257         enddo
 258       endif
 259 c
 260 c Compare angles
 261 c
 262       angn=0.0d0
 263       anorm=0
 264       nn = 0
 265       fract = 1.0d0
 266       npart = npiece_c
 267       do i=1,npart
 268         ishifc=ishift_c(i)
 269         nbeg = ifrag_c(1,i) + 3 - ishifc
 270         if (nbeg.lt.nn4) nbeg=nn4
 271         nend = ifrag_c(2,i)  - ishifc + 1
 272         if (nend.gt.nne) nend=nne
 273         if (nend.ge.nbeg) then
 274         nn = nn + nend - nbeg + 1
 275         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
 276      &    " nn",nn," ishifc",ishifc
 277         if (lprn) write (iout,*) "angles"
 278         longest=0
 279         ll = 0
 280         do j=nbeg,nend
 281 c          deltang = pinorm(phi(j)-phi_ref(j+ishifc,iprot))
 282           deltang=spherang(phi_ref(j+ishifc,iprot),
 283      &      theta_ref(j-1+ishifc,iprot),theta_ref(j+ishifc,iprot),
 284      &      phi(j),theta(j-1),theta(j))
 285           if (dabs(deltang).gt.diffang_max) then
 286             if (ll.gt.longest) longest = ll
 287             ll = 0
 288           else
 289             ll=ll+1
 290           endif
 291           if (ll.gt.longest) longest = ll
 292           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
 293      &     rad2deg*phi_ref(j+ishifc,iprot),rad2deg*deltang
 294           angn=angn+dabs(deltang)
 295         enddo
 296         longest=longest+3
 297         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
 298         if (lprn) write (iout,*)"segment",i," longest fragment within",
 299      &    diffang_max*rad2deg,":",longest," fraction",ff
 300         if (ff.lt.fract) fract = ff
 301         endif
 302       enddo
 303       if (nn.gt.0) anorm = angn/nn
 304       if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
 305       return
 306       end
 307 c-------------------------------------------------------------------------
 308       double precision function angnorm1(nang_pair,iang_pair,iprot,lprn)
 309       implicit none
 310       include 'DIMENSIONS'
 311       include 'DIMENSIONS.ZSCOPT'
 312       include 'DIMENSIONS.COMPAR'
 313       include 'COMMON.IOUNITS'
 314       include 'COMMON.VAR'
 315       include 'COMMON.COMPAR'
 316       include 'COMMON.CHAIN'
 317       include 'COMMON.GEO'
 318       logical lprn
 319       integer nang_pair,iang_pair(2,maxres)
 320       integer j,ia1,ia2
 321       double precision angn,deltang
 322       integer iprot
 323       double precision pinorm,spherang
 324       angn=0.0d0
 325       if (lprn) write (iout,'(80(1h*))')
 326       if (lprn) write (iout,*) "nang_pair",nang_pair
 327       if (lprn) write (iout,*) "angles"
 328       do j=1,nang_pair
 329         ia1 = iang_pair(1,j)
 330         ia2 = iang_pair(2,j)
 331 c        deltang = pinorm(phi(ia1)-phi_ref(ia2,iprot))
 332         deltang=spherang(phi_ref(ia2,iprot),theta_ref(ia2-1,iprot),
 333      &      theta_ref(ia2,iprot),phi(ia2),theta(ia2-1),theta(ia2))
 334         if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),
 335      &   rad2deg*phi_ref(ia2,iprot),rad2deg*deltang
 336         angn=angn+dabs(deltang)
 337       enddo
 338       if (lprn)
 339      &write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
 340       angnorm1 = angn/nang_pair
 341       return
 342       end
 343 c------------------------------------------------------------------------------
 344       subroutine angnorm12(diff,iprot)
 345       implicit none
 346       include 'DIMENSIONS'
 347       include 'DIMENSIONS.ZSCOPT'
 348       include 'DIMENSIONS.COMPAR'
 349       include 'COMMON.IOUNITS'
 350       include 'COMMON.VAR'
 351       include 'COMMON.COMPAR'
 352       include 'COMMON.CHAIN'
 353       include 'COMMON.GEO'
 354       integer iprot,nn4,nne,j
 355       double precision diff
 356       double precision pinorm,spherang
 357       diff=0.0d0
 358       nn4 = nstart_sup(iprot)+3
 359       nne = min0(nend_sup(iprot),nres)
 360 c      do j=nn4-1,nne
 361 c        diff = diff+dabs(pinorm(theta(j)-theta_ref(j,iprot)))
 362 c      enddo
 363       do j=nn4,nne
 364 c        diff = diff+dabs(pinorm(phi(j)-phi_ref(j,iprot)))
 365          diff=diff+spherang(phi_ref(j,iprot),theta_ref(j-1,iprot),
 366      &  theta_ref(j,iprot),phi(j),theta(j-1),theta(j))
 367       enddo
 368       diff=diff/(2*(nne-nn4)+3)
 369       return
 370       end
 371 c--------------------------------------------------------------------------------
 372       double precision function spherang(gam1,theta11,theta12,
 373      &   gam2,theta21,theta22)
 374       implicit none
 375       double precision gam1,theta11,theta12,gam2,theta21,theta22,
 376      &  x1,x2,xmed,f1,f2,fmed
 377       double precision tolx /1.0d-4/, tolf /1.0d-4/
 378       double precision sumcos
 379       double precision arcos,pinorm,sumangp
 380       integer it,maxit /100/
 381 c Calculate the difference of the angles of two superposed 4-redidue fragments
 382 c
 383 c       O      P
 384 c        \    /
 385 c     O'--C--C
 386 c             \
 387 c              P'
 388 c
 389 c The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
 390 c to achieve the minimum difference between the O'-C-O and P-C-P angles;
 391 c the sum of these angles is the difference returned by the function.
 392 c
 393 c 4/28/04 AL
 394 c If thetas match, take the difference of gamma and exit.
 395       if (dabs(theta11-theta12).lt.tolx
 396      & .and. dabs(theta21-theta22).lt.tolx) then
 397         spherang=dabs(pinorm(gam2-gam1))
 398         return
 399       endif
 400 c If the gammas are the same, take the difference of thetas and exit.
 401       x1=0.0d0
 402       x2=0.5d0*pinorm(gam2-gam1)
 403       if (dabs(x2) .lt. tolx) then
 404         spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
 405         return
 406       else if (x2.lt.0.0d0) then
 407         x1=x2
 408         x2=0.0d0
 409       endif
 410 c Else apply regula falsi method to compute optimum overlap of the terminal Calphas
 411       f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
 412       f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
 413       do it=1,maxit
 414         xmed=x1-f1*(x2-x1)/(f2-f1)
 415         fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
 416 c        write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
 417         if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx)
 418      &       .and. dabs(fmed).lt.tolf ) then
 419           x1=xmed
 420           f1=fmed
 421           goto 10
 422         else if ( fmed*f1.lt.0.0d0 ) then
 423           x2=xmed
 424           f2=fmed
 425         else
 426           x1=xmed
 427           f1=fmed
 428         endif
 429       enddo
 430    10 continue
 431       spherang=arcos(dcos(theta11)*dcos(theta12)
 432      & +dsin(theta11)*dsin(theta12)*dcos(x1))+
 433      & arcos(dcos(theta21)*dcos(theta22)+
 434      & dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
 435       return
 436       end
 437 c--------------------------------------------------------------------------------
 438       double precision function sumangp(gam1,theta11,theta12,gam2,
 439      & theta21,theta22,fi)
 440       implicit none
 441       double precision gam1,theta11,theta12,gam2,theta21,theta22,fi,
 442      & cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,
 443      & cosd2
 444 c derivarive of the sum of the difference of the angles of a 4-residue fragment.
 445       double precision arcos
 446       cost11=dcos(theta11)
 447       cost12=dcos(theta12)
 448       cost21=dcos(theta21)
 449       cost22=dcos(theta22)
 450       sint11=dsin(theta11)
 451       sint12=dsin(theta12)
 452       sint21=dsin(theta21)
 453       sint22=dsin(theta22)
 454       cosd1=cost11*cost12+sint11*sint12*dcos(fi)
 455       cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
 456       sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1)
 457      & +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)
 458       return
 459       end