X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fangnorm.f;fp=source%2Fwham%2Fsrc-NEWSC%2Fangnorm.f;h=2d1794243af72229fe73def6c772fc84a170f7b0;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/wham/src-NEWSC/angnorm.f b/source/wham/src-NEWSC/angnorm.f new file mode 100755 index 0000000..2d17942 --- /dev/null +++ b/source/wham/src-NEWSC/angnorm.f @@ -0,0 +1,439 @@ + subroutine add_angpair(ici,icj,nang_pair,iang_pair) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + integer ici,icj,nang_pair,iang_pair(2,maxres) + integer i,ian1,ian2 +c write (iout,*) "add_angpair: ici",ici," icj",icj, +c & " nang_pair",nang_pair + ian1=ici+2 + if (ian1.lt.4 .or. ian1.gt.nres) return + ian2=icj+2 +c write (iout,*) "ian1",ian1," ian2",ian2 + if (ian2.lt.4 .or. ian2.gt.nres) return + do i=1,nang_pair + if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return + enddo + nang_pair=nang_pair+1 + iang_pair(1,nang_pair)=ian1 + iang_pair(2,nang_pair)=ian2 + return + end +c------------------------------------------------------------------------- + subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract, + & lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.COMPAR' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + double precision pinorm,deltang + logical lprn + if (lprn) write (iout,'(80(1h*))') + angn=0.0d0 + nn = 0 + fract = 1.0d0 + npart = npiece(jfrag,1) + nn4 = nstart_sup+3 + nne = min0(nend_sup,nres) + if (lprn) write (iout,*) "nn4",nn4," nne",nne + do i=1,npart + nbeg = ifrag(1,i,jfrag) + 3 - ishif1 + if (nbeg.lt.nn4) nbeg=nn4 + nend = ifrag(2,i,jfrag) + 1 - ishif2 + if (nend.gt.nne) nend=nne + if (nend.ge.nbeg) then + nn = nn + nend - nbeg + 1 + if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend, + & " nn",nn," ishift1",ishif1," ishift2",ishif2 + if (lprn) write (iout,*) "angles" + longest=0 + ll = 0 + do j=nbeg,nend +c deltang = pinorm(phi(j)-phi_ref(j+ishif1)) + deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1), + & theta_ref(j+ishif1),phi(j),theta(j-1),theta(j)) + if (dabs(deltang).gt.diffang_max) then + if (ll.gt.longest) longest = ll + ll = 0 + else + ll=ll+1 + endif + if (ll.gt.longest) longest = ll + if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j), + & rad2deg*phi_ref(j+ishif1),rad2deg*deltang + angn=angn+dabs(deltang) + enddo + longest=longest+3 + ff = dfloat(longest)/dfloat(nend - nbeg + 4) + if (lprn) write (iout,*)"segment",i," longest fragment within", + & diffang_max*rad2deg,":",longest," fraction",ff + if (ff.lt.fract) fract = ff + endif + enddo + if (nn.gt.0) then + angn = angn/nn + else + angn = dwapi + endif + if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn, + & " fract",fract + return + end +c------------------------------------------------------------------------- + subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn, + & diffang_max,anorm,fract) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.COMPAR' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + integer ncont,icont(2,ncont),longest + double precision anorm,diffang_max,fract + integer npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece) + double precision pinorm + logical lprn + if (lprn) write (iout,'(80(1h*))') +c +c Determine the segments for which angles will be compared +c + nn4 = nstart_sup+3 + nne = min0(nend_sup,nres) + if (lprn) write (iout,*) "nn4",nn4," nne",nne + npart=npiece(jfrag,1) + npiece_c=0 + do i=1,npart +c write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) + if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. + & icont(1,1).gt.ifrag(2,i,jfrag)) goto 11 + jstart=1 + do while (jstart.lt.ncont .and. + & icont(1,jstart).lt.ifrag(1,i,jfrag)) +c write (iout,*) "jstart",jstart," icont",icont(1,jstart), +c & " ifrag",ifrag(1,i,jfrag) + jstart=jstart+1 + enddo +c write (iout,*) "jstart",jstart," icont",icont(1,jstart), +c & " ifrag",ifrag(1,i,jfrag) + if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11 + npiece_c=npiece_c+1 + ic1=icont(1,jstart) + ifrag_c(1,npiece_c)=icont(1,jstart) + jend=ncont + do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag)) +c write (iout,*) "jend",jend," icont",icont(1,jend), +c & " ifrag",ifrag(2,i,jfrag) + jend=jend-1 + enddo +c write (iout,*) "jend",jend," icont",icont(1,jend), +c & " ifrag",ifrag(2,i,jfrag) + ic2=icont(1,jend) + ifrag_c(2,npiece_c)=icont(1,jend)+1 + ishift_c(npiece_c)=ishif1 +c write (iout,*) "1: i",i," jstart:",jstart," jend",jend, +c & " ic1",ic1," ic2",ic2, +c & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) + 11 continue + if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then + idi=1 + else + idi=-1 + endif +c write (iout,*) "idi",idi + if (idi.eq.1) then + if (icont(2,1).gt.ifrag(2,i,jfrag) .or. + & icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12 + jstart=1 + do while (jstart.lt.ncont .and. + & icont(2,jstart).lt.ifrag(1,i,jfrag)) +c write (iout,*) "jstart",jstart," icont",icont(2,jstart), +c & " ifrag",ifrag(1,i,jfrag) + jstart=jstart+1 + enddo +c write (iout,*) "jstart",jstart," icont",icont(2,jstart), +c & " ifrag",ifrag(1,i,jfrag) + if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12 + npiece_c=npiece_c+1 + ic1=icont(2,jstart) + ifrag_c(2,npiece_c)=icont(2,jstart)+1 + jend=ncont + do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag)) +c write (iout,*) "jend",jend," icont",icont(2,jend), +c & " ifrag",ifrag(2,i,jfrag) + jend=jend-1 + enddo +c write (iout,*) "jend",jend," icont",icont(2,jend), +c & " ifrag",ifrag(2,i,jfrag) + else if (idi.eq.-1) then + if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or. + & icont(2,1).lt.ifrag(1,i,jfrag)) goto 12 + jstart=ncont + do while (jstart.gt.ncont .and. + & icont(2,jstart).lt.ifrag(1,i,jfrag)) +c write (iout,*) "jstart",jstart," icont",icont(2,jstart), +c & " ifrag",ifrag(1,i,jfrag) + jstart=jstart-1 + enddo +c write (iout,*) "jstart",jstart," icont",icont(2,jstart), +c & " ifrag",ifrag(1,i,jfrag) + if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12 + npiece_c=npiece_c+1 + ic1=icont(2,jstart) + ifrag_c(2,npiece_c)=icont(2,jstart)+1 + jend=1 + do while (jend.lt.ncont .and. + & icont(2,jend).gt.ifrag(2,i,jfrag)) +c write (iout,*) "jend",jend," icont",icont(2,jend), +c & " ifrag",ifrag(2,i,jfrag) + jend=jend+1 + enddo +c write (iout,*) "jend",jend," icont",icont(2,jend), +c & " ifrag",ifrag(2,i,jfrag) + endif + ic2=icont(2,jend) + if (ic2.lt.ic1) then + iic = ic1 + ic1 = ic2 + ic2 = iic + endif +c write (iout,*) "2: i",i," ic1",ic1," ic2",ic2, +c & " jstart:",jstart," jend",jend, +c & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) + ifrag_c(1,npiece_c)=ic1 + ifrag_c(2,npiece_c)=ic2+1 + ishift_c(npiece_c)=ishif2 + 12 continue + enddo + if (lprn) then + write (iout,*) "Before merge: npiece_c",npiece_c + do i=1,npiece_c + write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i) + enddo + endif +c +c Merge overlapping segments (e.g., avoid splitting helices) +c + i=1 + do while (i .lt. npiece_c) + if (ishift_c(i).eq.ishift_c(i+1) .and. + & ifrag_c(2,i).gt.ifrag_c(1,i+1)) then + ifrag_c(2,i)=ifrag_c(2,i+1) + do j=i+1,npiece_c + ishift_c(j)=ishift_c(j+1) + ifrag_c(1,j)=ifrag_c(1,j+1) + ifrag_c(2,j)=ifrag_c(2,j+1) + enddo + npiece_c=npiece_c-1 + else + i=i+1 + endif + enddo + if (lprn) then + write (iout,*) "After merge: npiece_c",npiece_c + do i=1,npiece_c + write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i) + enddo + endif +c +c Compare angles +c + angn=0.0d0 + anorm=0 + nn = 0 + fract = 1.0d0 + npart = npiece_c + do i=1,npart + ishifc=ishift_c(i) + nbeg = ifrag_c(1,i) + 3 - ishifc + if (nbeg.lt.nn4) nbeg=nn4 + nend = ifrag_c(2,i) - ishifc + 1 + if (nend.gt.nne) nend=nne + if (nend.ge.nbeg) then + nn = nn + nend - nbeg + 1 + if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend, + & " nn",nn," ishifc",ishifc + if (lprn) write (iout,*) "angles" + longest=0 + ll = 0 + do j=nbeg,nend +c deltang = pinorm(phi(j)-phi_ref(j+ishifc)) + deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc), + & theta_ref(j+ishifc),phi(j),theta(j-1),theta(j)) + if (dabs(deltang).gt.diffang_max) then + if (ll.gt.longest) longest = ll + ll = 0 + else + ll=ll+1 + endif + if (ll.gt.longest) longest = ll + if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j), + & rad2deg*phi_ref(j+ishifc),rad2deg*deltang + angn=angn+dabs(deltang) + enddo + longest=longest+3 + ff = dfloat(longest)/dfloat(nend - nbeg + 4) + if (lprn) write (iout,*)"segment",i," longest fragment within", + & diffang_max*rad2deg,":",longest," fraction",ff + if (ff.lt.fract) fract = ff + endif + enddo + if (nn.gt.0) anorm = angn/nn + if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract + return + end +c------------------------------------------------------------------------- + double precision function angnorm1(nang_pair,iang_pair,lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.COMPAR' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + logical lprn + integer nang_pair,iang_pair(2,maxres) + double precision pinorm + angn=0.0d0 + if (lprn) write (iout,'(80(1h*))') + if (lprn) write (iout,*) "nang_pair",nang_pair + if (lprn) write (iout,*) "angles" + do j=1,nang_pair + ia1 = iang_pair(1,j) + ia2 = iang_pair(2,j) +c deltang = pinorm(phi(ia1)-phi_ref(ia2)) + deltang=spherang(phi_ref(ia2),theta_ref(ia2-1), + & theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2)) + if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1), + & rad2deg*phi_ref(ia2),rad2deg*deltang + angn=angn+dabs(deltang) + enddo + if (lprn) + &write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair + angnorm1 = angn/nang_pair + return + end +c------------------------------------------------------------------------------ + subroutine angnorm12(diff) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.COMPAR' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + double precision pinorm + diff=0.0d0 + nn4 = nstart_sup+3 + nne = min0(nend_sup,nres) +c do j=nn4-1,nne +c diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j))) +c enddo + do j=nn4,nne +c diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j))) + diff=diff+spherang(phi_ref(j),theta_ref(j-1), + & theta_ref(j),phi(j),theta(j-1),theta(j)) + enddo + return + end +c-------------------------------------------------------------------------------- + double precision function spherang(gam1,theta11,theta12, + & gam2,theta21,theta22) + implicit none + double precision gam1,theta11,theta12,gam2,theta21,theta22, + & x1,x2,xmed,f1,f2,fmed + double precision tolx /1.0d-4/, tolf /1.0d-4/ + double precision sumcos + double precision arcos,pinorm,sumangp + integer it,maxit /100/ +c Calculate the difference of the angles of two superposed 4-redidue fragments +c +c O P +c \ / +c O'--C--C +c \ +c P' +c +c The fragment O'-C-C-P' is rotated by angle fi about the C-C axis +c to achieve the minimum difference between the O'-C-O and P-C-P angles; +c the sum of these angles is the difference returned by the function. +c +c 4/28/04 AL +c If thetas match, take the difference of gamma and exit. + if (dabs(theta11-theta12).lt.tolx + & .and. dabs(theta21-theta22).lt.tolx) then + spherang=dabs(pinorm(gam2-gam1)) + return + endif +c If the gammas are the same, take the difference of thetas and exit. + x1=0.0d0 + x2=0.5d0*pinorm(gam2-gam1) + if (dabs(x2) .lt. tolx) then + spherang=dabs(theta11-theta21)+dabs(theta12-theta22) + return + else if (x2.lt.0.0d0) then + x1=x2 + x2=0.0d0 + endif +c Else apply regula falsi method to compute optimum overlap of the terminal Calphas + f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1) + f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2) + do it=1,maxit + xmed=x1-f1*(x2-x1)/(f2-f1) + fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed) +c write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed + if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx) + & .and. dabs(fmed).lt.tolf ) then + x1=xmed + f1=fmed + goto 10 + else if ( fmed*f1.lt.0.0d0 ) then + x2=xmed + f2=fmed + else + x1=xmed + f1=fmed + endif + enddo + 10 continue + spherang=arcos(dcos(theta11)*dcos(theta12) + & +dsin(theta11)*dsin(theta12)*dcos(x1))+ + & arcos(dcos(theta21)*dcos(theta22)+ + & dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1)) + return + end +c-------------------------------------------------------------------------------- + double precision function sumangp(gam1,theta11,theta12,gam2, + & theta21,theta22,fi) + implicit none + double precision gam1,theta11,theta12,gam2,theta21,theta22,fi, + & cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1, + & cosd2 +c derivarive of the sum of the difference of the angles of a 4-residue fragment. + double precision arcos + cost11=dcos(theta11) + cost12=dcos(theta12) + cost21=dcos(theta21) + cost22=dcos(theta22) + sint11=dsin(theta11) + sint12=dsin(theta12) + sint21=dsin(theta21) + sint22=dsin(theta22) + cosd1=cost11*cost12+sint11*sint12*dcos(fi) + cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi) + sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1) + & +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2) + return + end