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