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 #define DEBUG
  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.nz_start .or. ichain1.gt.nz_end .or.
  33      &          ichain2.lt.nz_start .or. ichain2.gt.nz_end) 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.nz_start .or. ichain1.gt.nz_end .or.
  58      &          ichain2.lt.nz_start .or. ichain2.gt.nz_end) 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       endif
  94 #ifdef DEBUG
  95       write (iout,*) "iperm",iperm," rms",rms
  96 #endif
  97       ENDDO
  98       rmscalc=rmsmin
  99 #ifdef DEBUG
 100       write (iout,*) "ipermmin",ipermmin," rmsmin",rmsmin
 101 #endif
 102       return
 103       end
 104 c------------------------------------------------------------------------
 105       double precision function rmscalc_thet(ttheta,theta_reff,
 106      &     iperm)
 107       implicit none
 108       include 'DIMENSIONS'
 109       include 'COMMON.IOUNITS'
 110       include 'COMMON.CHAIN'
 111       include 'COMMON.INTERACT'
 112       include 'COMMON.VAR'
 113       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 114       double precision ttheta(maxres),theta_reff(maxres),rmsthet,dtheta
 115       rmsthet = 0.0d0
 116       nnnn=0
 117       do ichain=1,nchain
 118         indchain=tabpermchain(ichain,iperm)
 119 c        write (iout,*) "ichain",ichain," iperm",iperm,
 120 c     &    " indchain",indchain
 121 c        call flush(iout)
 122         do k=3,chain_length(ichain)
 123           kchain1=chain_border(1,ichain)+k-1
 124           kchain2=chain_border(1,indchain)+k-1
 125           nnnn=nnnn+1
 126           dtheta = ttheta(kchain2)-theta_reff(kchain1)
 127 c                write (iout,*) k,theta(k),theta_ref(k,iref,ib,iprot),
 128 c     &            dtheta
 129           rmsthet = rmsthet+dtheta*dtheta
 130         enddo
 131       enddo
 132       nnnn=nnnn-1
 133       rmsthet=dsqrt(rmsthet/nnnn)
 134 #ifdef DEBUG
 135       write (iout,*) "nnnn",nnnn," rmsthet",rmsthet
 136 #endif
 137       rmscalc_thet=rmsthet
 138       return
 139       end
 140 c------------------------------------------------------------------------
 141       double precision function rmscalc_phi(pphi,phi_reff,iperm)
 142       implicit none
 143       include 'DIMENSIONS'
 144       include 'COMMON.IOUNITS'
 145       include 'COMMON.CHAIN'
 146       include 'COMMON.INTERACT'
 147       include 'COMMON.VAR'
 148       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 149       double precision pphi(maxres),phi_reff(maxres),rmsphi,dphi
 150       double precision pinorm
 151       rmsphi = 0.0d0
 152       nnnn=0
 153       do ichain=1,nchain
 154         indchain=tabpermchain(ichain,iperm)
 155         do k=4,chain_length(ichain)
 156           kchain1=chain_border(1,ichain)+k-1
 157           kchain2=chain_border(1,indchain)+k-1
 158           nnnn=nnnn+1
 159           dphi=pinorm(pphi(kchain2)-phi_reff(kchain1))
 160 c         write (iout,*) k,phi(k),phi_ref(k,iref,ib,iprot),
 161 c     &   pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
 162           rmsphi = rmsphi + dphi*dphi
 163         enddo
 164       enddo
 165       nnnn=nnnn-1
 166       rmsphi=dsqrt(rmsphi/nnnn)
 167 #ifdef DEBUG
 168       write (iout,*) "nnnn",nnnn," rmsphi",rmsphi
 169 #endif
 170       rmscalc_phi=rmsphi
 171       return
 172       end
 173 c------------------------------------------------------------------------
 174       double precision function rmscalc_side(xxtabb,yytabb,zztabb,
 175      & xxreff,yyreff,zzreff,iperm)
 176       implicit none
 177       include 'DIMENSIONS'
 178       include 'COMMON.IOUNITS'
 179       include 'COMMON.CHAIN'
 180       include 'COMMON.INTERACT'
 181       include 'COMMON.VAR'
 182       integer iperm,k,ichain,indchain,kchain1,kchain2,nnnn
 183       double precision xxtabb(maxres),yytabb(maxres),zztabb(maxres),
 184      & xxreff(maxres),yyreff(maxres),zzreff(maxres),rmsside,
 185      & dxref,dyref,dzref
 186       rmsside = 0.0d0
 187       nnnn=0
 188       do ichain=1,nchain
 189         indchain=tabpermchain(ichain,iperm)
 190         do k=1,chain_length(ichain)
 191           kchain1=chain_border(1,ichain)+k-1
 192           kchain2=chain_border(1,indchain)+k-1
 193           if (itype(kchain1).eq.ntyp1) cycle
 194           nnnn=nnnn+1
 195           dxref = xxtabb(kchain2)-xxreff(kchain1)
 196           dyref = yytabb(kchain2)-yyreff(kchain1)
 197           dzref = zztabb(kchain2)-zzreff(kchain1)
 198           rmsside = rmsside + dxref*dxref+dyref*dyref+dzref*dzref
 199         enddo
 200       enddo
 201       rmsside=dsqrt(rmsside/nnnn)
 202 #ifdef DEBUG
 203       write (iout,*) "nnnn",nnnn," rmsside",rmsside
 204 #endif
 205 #undef DEBUG
 206       rmscalc_side=rmsside
 207       return
 208       end