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