source/maxlik/src-Fmatch/rmscalc.F

   1       double precision function rmsnat(ib,jcon,iref,iprot,iperm)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'COMMON.IOUNITS'
   6       include 'COMMON.COMPAR'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.VAR'
  10       double precision creff(3,maxres2),cc(3,maxres2)
  11       integer jcon,ib,iref,iprot,iperm
  12       double precision rmscalc
  13       rmsnat=rmscalc(c(1,1),cref(1,1,iref,ib,iprot),jcon,iref,iprot,
  14      &  iperm)
  15       return
  16       end
  17 c------------------------------------------------------------------------
  18       double precision function rmscalc(ccc,cccref,jcon,iref,iprot,
  19      &   ipermmin)
  20       implicit none
  21       include 'DIMENSIONS'
  22       include 'DIMENSIONS.ZSCOPT'
  23       include 'COMMON.IOUNITS'
  24       include 'COMMON.COMPAR'
  25       include 'COMMON.CHAIN'
  26       include 'COMMON.INTERACT'
  27       include 'COMMON.VAR'
  28       double precision cccref(3,maxres2),creff(3,maxres2),
  29      &  ccc(3,maxres2),cc(3,maxres2)
  30       double precision przes(3),obrot(3,3)
  31       logical non_conv
  32       integer i,ii,j,jcon,ib,iref,iprot,ichain,indchain,ichain1,ichain2,
  33      &  iperm,ipermmin
  34       double precision rms,rmsmin
  35 C Loop over chain permutations
  36       rmsmin=1.0d10
  37       DO IPERM=1,NPERMCHAIN
  38       ii=0
  39       if (.not.sconly(iprot)) then
  40         do ichain=1,nchain
  41           indchain=tabpermchain(ichain,iperm)
  42 #ifdef DEBUG
  43           write (iout,*) "ichain",ichain," indchain",indchain
  44           write (iout,*) "chain_border",chain_border(1,ichain),
  45      &      chain_border(2,ichain)
  46 #endif
  47           do i=1,chain_length(ichain)
  48             ichain1=chain_border(1,ichain)+i-1
  49             ichain2=chain_border(1,indchain)+i-1
  50             ii=ii+1
  51 #ifdef DEBUG
  52             write (iout,*) "back",ii," ichain1",ichain1,
  53      &       " ichain2",ichain2," i",i,chain_border(1,ichain)+i-1
  54 #endif
  55             do j=1,3
  56               cc(j,ii)=ccc(j,ichain2)
  57               creff(j,ii)=cccref(j,ichain1)
  58             enddo
  59 #ifdef DEBUG
  60             write (iout,'(3f10.5,5x,3f10.5)')
  61      &       (cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
  62 #endif
  63           enddo
  64         enddo
  65       endif
  66       if (.not.caonly(iprot)) then
  67         do ichain=1,nchain
  68           indchain=tabpermchain(ichain,iperm)
  69           do i=1,chain_length(ichain)
  70             ichain1=chain_border(1,ichain)+i-1
  71             ichain2=chain_border(1,indchain)+i-1
  72             if (itype(ichain1).ne.10) then
  73               ii=ii+1
  74 #ifdef DEBUG
  75               write (iout,*) "side",ii," ichain1",ichain1,
  76      &         " ichain2",ichain2
  77 #endif
  78               do j=1,3
  79                 cc(j,ii)=ccc(j,ichain2+nres)
  80                 creff(j,ii)=cccref(j,ichain1+nres)
  81               enddo
  82 #ifdef DEBUG
  83               write (iout,'(3f10.5,5x,3f10.5)')
  84      &        (cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
  85 #endif
  86             endif
  87           enddo
  88         enddo
  89       endif
  90 c      write (iout,*) "rmscalc: iprot",iprot," nsup",nsup(iprot)," ii",ii
  91       call fitsq(rms,cc(1,1),creff(1,1),ii,przes,obrot,non_conv)
  92       if (non_conv) then
  93         write (iout,*) 'Error: FITSQ non-convergent, iprot',iprot,
  94      &   ' jcon',jcon,' iref',iref
  95         rms=1.0d2
  96       else if (rms.lt.-1.0d-6) then
  97         print *,'Error: rms^2 = ',rms,jcon
  98         rms = 1.0d2
  99       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
 100         rmscalc=0.0d0
 101       else
 102         rms = dsqrt(rms)
 103       endif
 104       if (rms.lt.rmsmin) then
 105         rmsmin=rms
 106         ipermmin=iperm
 107       endif
 108 #ifdef DEBUG
 109       write (iout,*) "iperm",iperm," rms",rms
 110 #endif
 111       ENDDO
 112       rmscalc=rmsmin
 113 #ifdef DEBUG
 114       write (iout,*) "ipermmin",ipermmin," rmsmin",rmsmin
 115 #endif
 116       return
 117       end
 118 c------------------------------------------------------------------------
 119       double precision function rmscalc_thet(ttheta,theta_reff,
 120      &     iperm)
 121       implicit none
 122       include 'DIMENSIONS'
 123       include 'DIMENSIONS.ZSCOPT'
 124       include 'COMMON.IOUNITS'
 125       include 'COMMON.COMPAR'
 126       include 'COMMON.CHAIN'
 127       include 'COMMON.INTERACT'
 128       include 'COMMON.VAR'
 129       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 130       double precision ttheta(maxres),theta_reff(maxres),rmsthet,dtheta
 131       rmsthet = 0.0d0
 132       nnnn=0
 133       do ichain=1,nchain
 134         indchain=tabpermchain(ichain,iperm)
 135 c        write (iout,*) "ichain",ichain," iperm",iperm,
 136 c     &    " indchain",indchain
 137         call flush(iout)
 138         do k=3,chain_length(ichain)
 139           kchain1=chain_border(1,ichain)+k-1
 140           kchain2=chain_border(1,indchain)+k-1
 141           nnnn=nnnn+1
 142           dtheta = ttheta(kchain2)-theta_reff(kchain1)
 143 c                write (iout,*) k,theta(k),theta_ref(k,iref,ib,iprot),
 144 c     &            dtheta
 145           rmsthet = rmsthet+dtheta*dtheta
 146         enddo
 147       enddo
 148       nnnn=nnnn-1
 149       rmsthet=dsqrt(rmsthet/nnnn)
 150 #ifdef DEBUG
 151       write (iout,*) "nnnn",nnnn," rmsthet",rmsthet
 152 #endif
 153       rmscalc_thet=rmsthet
 154       return
 155       end
 156 c------------------------------------------------------------------------
 157       double precision function rmscalc_phi(pphi,phi_reff,iperm)
 158       implicit none
 159       include 'DIMENSIONS'
 160       include 'DIMENSIONS.ZSCOPT'
 161       include 'COMMON.IOUNITS'
 162       include 'COMMON.COMPAR'
 163       include 'COMMON.CHAIN'
 164       include 'COMMON.INTERACT'
 165       include 'COMMON.VAR'
 166       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 167       double precision pphi(maxres),phi_reff(maxres),rmsphi,dphi
 168       double precision pinorm
 169       rmsphi = 0.0d0
 170       nnnn=0
 171       do ichain=1,nchain
 172         indchain=tabpermchain(ichain,iperm)
 173         do k=4,chain_length(ichain)
 174           kchain1=chain_border(1,ichain)+k-1
 175           kchain2=chain_border(1,indchain)+k-1
 176           nnnn=nnnn+1
 177           dphi=pinorm(pphi(kchain2)-phi_reff(kchain1))
 178 c         write (iout,*) k,phi(k),phi_ref(k,iref,ib,iprot),
 179 c     &   pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
 180           rmsphi = rmsphi + dphi*dphi
 181         enddo
 182       enddo
 183       nnnn=nnnn-1
 184       rmsphi=dsqrt(rmsphi/nnnn)
 185 #ifdef DEBUG
 186       write (iout,*) "nnnn",nnnn," rmsphi",rmsphi
 187 #endif
 188       rmscalc_phi=rmsphi
 189       return
 190       end
 191 c------------------------------------------------------------------------
 192       double precision function rmscalc_side(xxtabb,yytabb,zztabb,
 193      & xxreff,yyreff,zzreff,iperm)
 194       implicit none
 195       include 'DIMENSIONS'
 196       include 'DIMENSIONS.ZSCOPT'
 197       include 'COMMON.IOUNITS'
 198       include 'COMMON.COMPAR'
 199       include 'COMMON.CHAIN'
 200       include 'COMMON.INTERACT'
 201       include 'COMMON.VAR'
 202       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 203       double precision xxtabb(maxres),yytabb(maxres),zztabb(maxres),
 204      & xxreff(maxres),yyreff(maxres),zzreff(maxres),rmsside,
 205      & dxref,dyref,dzref
 206       rmsside = 0.0d0
 207       nnnn=0
 208       do ichain=1,nchain
 209         indchain=tabpermchain(ichain,iperm)
 210         do k=1,chain_length(ichain)
 211           kchain1=chain_border(1,ichain)+k-1
 212           kchain2=chain_border(1,indchain)+k-1
 213           if (itype(kchain1).eq.ntyp1) cycle
 214           nnnn=nnnn+1
 215           dxref = xxtabb(kchain2)-xxreff(kchain1)
 216           dyref = yytabb(kchain2)-yyreff(kchain1)
 217           dzref = zztabb(kchain2)-zzreff(kchain1)
 218           rmsside = rmsside + dxref*dxref+dyref*dyref+dzref*dzref
 219         enddo
 220       enddo
 221       rmsside=dsqrt(rmsside/nnnn)
 222 #ifdef DEBUG
 223       write (iout,*) iii,iref," nnnn",nnnn," rmsside",rmsside
 224 #endif
 225       rmscalc_side=rmsside
 226       return
 227       end