source/unres/src-HCD-5D/rmscalc.F

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