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