source/cluster/wham/src-M-SAXS-homology/rmscalc.F

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