source/wham/src-M-SAXS-homology/angnorm.f

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