source/wham/src-M/rmscalc.f

   1       double precision function rmscalc(ishif,i,j,jcon,lprn)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.COMPAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.COMPAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CONTROL'
  12       double precision przes(3),obrot(3,3)
  13       double precision creff(3,maxres2),cc(3,maxres2)
  14       logical iadded(maxres)
  15       integer inumber(2,maxres)
  16       common /ccc/ creff,cc,iadded,inumber
  17       logical lprn
  18       logical non_conv
  19       integer ishif,i,j
  20       if (lprn) then
  21         write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
  22         write (iout,*) "npiece",npiece(j,i)
  23         call flush(iout)
  24       endif
  25 c      write (iout,*) "symetr",symetr
  26 c      call flush(iout)
  27       nperm=1
  28       do idup=1,symetr
  29       nperm=nperm*idup
  30       enddo
  31 c      write (iout,*) "nperm",nperm
  32 c      call flush(iout)
  33       do kkk=1,nperm
  34       idup=0
  35       do l=1,nres
  36         iadded(l)=.false.
  37       enddo
  38 c      write (iout,*) "kkk",kkk
  39 c      call flush(iout)
  40       do k=1,npiece(j,i)
  41         if (i.eq.1) then
  42           if (lprn) then
  43             write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",
  44      &         ifrag(1,k,j),ifrag(2,k,j)
  45             call flush(iout)
  46           endif
  47           call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
  48 c          write (iout,*) "Exit cprep"
  49 c          call flush(iout)
  50 c          write (iout,*) "ii=",ii
  51         else
  52           kk = ipiece(k,j,i)
  53 c          write (iout,*) "kk",kk," npiece",npiece(kk,1)
  54           do l=1,npiece(kk,1)
  55             if (lprn) then
  56               write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,
  57      &          " l=",l," adding fragment",
  58      &          ifrag(1,l,kk),ifrag(2,l,kk)
  59               call flush(iout)
  60             endif
  61             call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
  62 c            write (iout,*) "After cprep"
  63 c            call flush(iout)
  64           enddo
  65         endif
  66       enddo
  67       enddo
  68       if (lprn) then
  69         write (iout,*) "tuszukaj"
  70         do kkk=1,nperm
  71           do k=1,idup
  72             write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),
  73      &        inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
  74           enddo
  75         enddo
  76         call flush(iout)
  77       endif
  78       rminrms=1.0d10
  79       do kkk=1,nperm
  80       call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
  81       if (non_conv) then
  82         print *,'Error: FITSQ non-convergent, jcon',jcon,i
  83         rms = 1.0d10
  84       else if (rms.lt.-1.0d-6) then
  85         print *,'Error: rms^2 = ',rms,jcon,i
  86         rms = 1.0d10
  87       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
  88         rms = 0.0d0
  89       endif
  90 c      write (iout,*) "rmsmin", rminrms, "rms", rms
  91       if (rms.le.rminrms) rminrms=rms
  92       enddo
  93       rmscalc = dsqrt(rminrms)
  94 c      write (iout, *) "analysys", rmscalc,anatemp
  95       return
  96       end
  97 c-------------------------------------------------------------------------
  98       subroutine cprep(if1,if2,ishif,idup,kwa)
  99       implicit real*8 (a-h,o-z)
 100       include 'DIMENSIONS'
 101       include 'DIMENSIONS.ZSCOPT'
 102       include 'DIMENSIONS.COMPAR'
 103       include 'COMMON.CONTROL'
 104       include 'COMMON.IOUNITS'
 105       include 'COMMON.COMPAR'
 106       include 'COMMON.CHAIN'
 107       include 'COMMON.INTERACT'
 108       include 'COMMON.VAR'
 109       double precision przes(3),obrot(3,3)
 110       double precision creff(3,maxres2),cc(3,maxres2)
 111       logical iadded(maxres)
 112       integer inumber(2,maxres),iistrart,kwa,blar
 113       common /ccc/ creff,cc,iadded,inumber
 114 c      write (iout,*) "Calling cprep symetr",symetr," kwa",kwa
 115       nperm=1
 116       do blar=1,symetr
 117       nperm=nperm*blar
 118       enddo
 119 c      write (iout,*) "nperm",nperm
 120       kkk=kwa
 121 c      ii=0
 122       do l=if1,if2
 123 c        write (iout,*) "l",l," iadded",iadded(l)
 124 c        call flush(iout)
 125         if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l))
 126      &  then
 127           idup=idup+1
 128           iadded(l)=.true.
 129           inumber(1,idup)=l
 130           inumber(2,idup)=l+ishif
 131           do m=1,3
 132             creff(m,idup)=cref(m,l,kkk)
 133             cc(m,idup)=c(m,l+ishif)
 134           enddo
 135         endif
 136       enddo
 137       return
 138       end
 139 c-------------------------------------------------------------------------
 140       double precision function rmsnat(jcon)
 141       implicit real*8 (a-h,o-z)
 142       include 'DIMENSIONS'
 143       include 'DIMENSIONS.ZSCOPT'
 144       include 'DIMENSIONS.COMPAR'
 145       include 'COMMON.IOUNITS'
 146       include 'COMMON.COMPAR'
 147       include 'COMMON.CHAIN'
 148       include 'COMMON.INTERACT'
 149       include 'COMMON.VAR'
 150       include 'COMMON.CONTROL'
 151       double precision przes(3),obrot(3,3),cc(3,2*maxres),
 152      &   ccref(3,2*maxres)
 153       logical non_conv
 154       integer ishif,i,j,resprzesun
 155       rminrms=10.0d10
 156       rmsminsing=10d10
 157       nperm=1
 158 C      write (iout,*) "tu2", nres,nsup
 159       noverlap=nres
 160       if (nres.gt.nsup) noverlap=nsup
 161       write (iout,*) "tu3,",noverlap
 162       do i=1,symetr
 163        nperm=nperm*i
 164       enddo
 165       do kkk=1,nperm
 166        nnsup=0
 167        do i=1,noverlap
 168         if (itype(i).ne.ntyp1) then
 169           nnsup=nnsup+1
 170           do j=1,3
 171             cc(j,nnsup)=c(j,i)
 172             ccref(j,nnsup)=cref(j,i,kkk)
 173           enddo
 174         endif
 175        enddo
 176        call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv)
 177        if (non_conv) then
 178         print *,'Error: FITSQ non-convergent, jcon',jcon,i
 179         rms=1.0d10
 180        else if (rms.lt.-1.0d-6) then
 181         print *,'Error: rms^2 = ',rms,jcon,i
 182         rms = 1.0d10
 183        else if (rms.ge.1.0d-6 .and. rms.lt.0) then
 184         rms=0.0d0
 185        endif
 186        if (rms.le.rminrms) rminrms=rms
 187 c       write (iout,*) "kkk",kkk," rmsnat",rms , rminrms
 188       enddo
 189       rmsnat = dsqrt(rminrms)
 190 C      write (iout,*)  "analysys",rmsnat, anatemp
 191 C      liczenie rmsdla pojedynczego lancucha
 192       return
 193       end
 194 c-----------------------------------------------------------------------------
 195       double precision function gyrate(jcon)
 196       implicit real*8 (a-h,o-z)
 197       include 'DIMENSIONS'
 198       include 'COMMON.INTERACT'
 199       include 'COMMON.CHAIN'
 200       double precision cen(3),rg
 201
 202       do j=1,3
 203        cen(j)=0.0d0
 204       enddo
 205
 206       do i=nnt,nct
 207           do j=1,3
 208             cen(j)=cen(j)+c(j,i)
 209           enddo
 210       enddo
 211       do j=1,3
 212             cen(j)=cen(j)/dble(nct-nnt+1)
 213       enddo
 214       rg = 0.0d0
 215       do i = nnt, nct
 216         do j=1,3
 217          rg = rg + (c(j,i)-cen(j))**2
 218         enddo
 219       end do
 220       gyrate = dsqrt(rg/dble(nct-nnt+1))
 221       return
 222       end