source/wham/src-M-SAXS/rmscalc.F

   1       double precision function rmscalc_frag(ishif,i,j,jcon,kkk,
   2      &  lprn)
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'DIMENSIONS.ZSCOPT'
   6       include 'DIMENSIONS.COMPAR'
   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       integer kkk
  22       if (lprn) then
  23         write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
  24         write (iout,*) "npiece",npiece(j,i)
  25         write (iout,*)  "kkk",kkk
  26         call flush(iout)
  27       endif
  28 c      write (iout,*) "symetr",symetr
  29 c      call flush(iout)
  30 c      nperm=1
  31 c      do idup=1,symetr
  32 c      nperm=nperm*idup
  33 c      enddo
  34 c      write (iout,*) "nperm",nperm
  35 c      call flush(iout)
  36 c      do kkk=1,nperm
  37       idup=0
  38       do l=1,nres
  39         iadded(l)=.false.
  40       enddo
  41 c      write (iout,*) "kkk",kkk
  42 c      call flush(iout)
  43       do k=1,npiece(j,i)
  44         if (i.eq.1) then
  45           if (lprn) then
  46             write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",
  47      &         ifrag(1,k,j),ifrag(2,k,j)
  48             call flush(iout)
  49           endif
  50           call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
  51 c          write (iout,*) "Exit cprep"
  52 c          call flush(iout)
  53 c          write (iout,*) "ii=",ii
  54         else
  55           kk = ipiece(k,j,i)
  56 c          write (iout,*) "kk",kk," npiece",npiece(kk,1)
  57           do l=1,npiece(kk,1)
  58             if (lprn) then
  59               write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,
  60      &          " l=",l," adding fragment",
  61      &          ifrag(1,l,kk),ifrag(2,l,kk)
  62               call flush(iout)
  63             endif
  64             call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
  65 c            write (iout,*) "After cprep"
  66 c            call flush(iout)
  67           enddo
  68         endif
  69       enddo
  70       if (lprn) then
  71         write (iout,*) "tuszukaj"
  72 c        do kkk=1,nperm
  73           do k=1,idup
  74             write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),
  75      &        inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
  76           enddo
  77 c        enddo
  78         call flush(iout)
  79       endif
  80 c      rminrms=1.0d10
  81 c      do kkk=1,nperm
  82       call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
  83       if (non_conv) then
  84         print *,'Error: FITSQ non-convergent, jcon',jcon,i
  85         rms = 1.0d10
  86       else if (rms.lt.-1.0d-6) then
  87         print *,'Error: rms^2 = ',rms,jcon,i
  88         rms = 1.0d10
  89       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
  90         rms = 0.0d0
  91       endif
  92 c      write (iout,*) "rmsmin", rminrms, "rms", rms
  93 c      if (rms.le.rminrms) rminrms=rms
  94 c      enddo
  95       rmscalc_frag = dsqrt(rms)
  96 c      write (iout, *) "analysys", rmscalc,anatemp
  97       return
  98       end
  99 c-------------------------------------------------------------------------
 100       subroutine cprep(if1,if2,ishif,idup,kwa)
 101       implicit real*8 (a-h,o-z)
 102       include 'DIMENSIONS'
 103       include 'DIMENSIONS.ZSCOPT'
 104       include 'DIMENSIONS.COMPAR'
 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       integer ll,iperm
 117 c      write (iout,*) "Calling cprep if1",if1," if2",if2," ishif",ishif,
 118 c     &    " kwa",kwa
 119 c      nperm=1
 120 c      do blar=1,symetr
 121 c      nperm=nperm*blar
 122 c      enddo
 123 c      write (iout,*) "nperm",nperm
 124 c      kkk=kwa
 125 c      ii=0
 126       do l=if1,if2
 127 c        write (iout,*) "l",l," iadded",iadded(l)," ireschain",
 128 c     &     ireschain(l),ireschain(l+ishif)
 129 c        call flush(iout)
 130         if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)
 131      &   .and. ireschain(l+ishif).gt.0 .and. ireschain(l).gt.0 .and.
 132      &   ireschain(l).eq.ireschain(l+ishif)) then
 133           idup=idup+1
 134           iadded(l)=.true.
 135           inumber(1,idup)=l
 136           inumber(2,idup)=l+ishif
 137           ll=iperm(l+ishif,kwa)
 138           do m=1,3
 139             creff(m,idup)=cref(m,l)
 140             cc(m,idup)=c(m,ll)
 141           enddo
 142 c          write (iout,'(2i5,3f10.5,5x,3f10.5)') l,ll,
 143 c     &     (creff(m,idup),m=1,3),(cc(m,idup),m=1,3)
 144         endif
 145       enddo
 146 c      write (iout,*) "idup",idup
 147       return
 148       end
 149 c-------------------------------------------------------------------------
 150       double precision function rmsnat(jcon,ipermmin)
 151       implicit real*8 (a-h,o-z)
 152       include 'DIMENSIONS'
 153       include 'DIMENSIONS.ZSCOPT'
 154       include 'DIMENSIONS.COMPAR'
 155       include 'COMMON.IOUNITS'
 156       include 'COMMON.COMPAR'
 157       include 'COMMON.CHAIN'
 158       include 'COMMON.INTERACT'
 159       include 'COMMON.VAR'
 160       include 'COMMON.CONTROL'
 161       integer ipermmin
 162       rmsnat = rmscalc(c(1,1),cref(1,1),ipermmin)
 163       return
 164       end
 165 c-----------------------------------------------------------------------------
 166       double precision function rmscalc(ccc,cccref,ipermmin)
 167       implicit none
 168       include 'DIMENSIONS'
 169       include 'COMMON.IOUNITS'
 170       include 'COMMON.CHAIN'
 171       include 'COMMON.INTERACT'
 172       include 'COMMON.VAR'
 173       double precision cccref(3,maxres2),creff(3,maxres2),
 174      &  ccc(3,maxres2),cc(3,maxres2)
 175       double precision przes(3),obrot(3,3)
 176       logical non_conv
 177       integer i,ii,j,ib,ichain,indchain,ichain1,ichain2,
 178      &  iperm,ipermmin
 179       double precision rms,rmsmin
 180 C Loop over chain permutations
 181 c      write (iout,*) "iz_sc",iz_sc
 182       rmsmin=1.0d10
 183       DO IPERM=1,NPERMCHAIN
 184 c      write (iout,*) "iperm",iperm
 185       ii=0
 186       if (iz_sc.lt.2) then
 187         do ichain=1,nchain
 188           indchain=tabpermchain(ichain,iperm)
 189 #ifdef DEBUG
 190           write (iout,*) "ichain",ichain," indchain",indchain
 191           write (iout,*) "chain_border",chain_border(1,ichain),
 192      &      chain_border(2,ichain)
 193 #endif
 194           do i=1,chain_length(ichain)
 195 c          do i=nstart_sup(ichain),nend_sup(ichain)
 196             ichain1=chain_border(1,ichain)+i-1
 197             ichain2=chain_border(1,indchain)+i-1
 198             if (ichain1.lt.nstart_sup .or. ichain1.gt.nend_sup .or.
 199      &          ichain2.lt.nstart_sup .or. ichain2.gt.nend_sup) cycle
 200             ii=ii+1
 201 #ifdef DEBUG
 202             write (iout,*) "back",ii," ichain1",ichain1,
 203      &       " ichain2",ichain2," i",i,chain_border(1,ichain)+i-1
 204 #endif
 205             do j=1,3
 206               cc(j,ii)=ccc(j,ichain2)
 207               creff(j,ii)=cccref(j,ichain1)
 208             enddo
 209 #ifdef DEBUG
 210             write (iout,'(2i5,3f10.5,5x,3f10.5)')
 211      &       ichain1,ii,(cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
 212 #endif
 213           enddo
 214         enddo
 215       endif
 216       if (iz_sc.gt.0) then
 217         do ichain=1,nchain
 218           indchain=tabpermchain(ichain,iperm)
 219           do i=1,chain_length(ichain)
 220 c          do i=nstart_sup(ichain),nend_sup(ichain)
 221             ichain1=chain_border(1,ichain)+i-1
 222             ichain2=chain_border(1,indchain)+i-1
 223             if (ichain1.lt.nstart_sup .or. ichain1.gt.nend_sup .or.
 224      &          ichain2.lt.nstart_sup .or. ichain2.gt.nend_sup) cycle
 225             if (itype(ichain1).ne.10) then
 226               ii=ii+1
 227 #ifdef DEBUG
 228               write (iout,*) "side",ii," ichain1",ichain1,
 229      &         " ichain2",ichain2
 230 #endif
 231               do j=1,3
 232                 cc(j,ii)=ccc(j,ichain2+nres)
 233                 creff(j,ii)=cccref(j,ichain1+nres)
 234               enddo
 235 #ifdef DEBUG
 236               write (iout,'(2i5,3f10.5,5x,3f10.5)')
 237      &        ichain1+nres,ii,(cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
 238 #endif
 239             endif
 240           enddo
 241         enddo
 242       endif
 243 c      write (iout,*) "rmscalc: iprot",iprot," nsup",nsup(iprot)," ii",ii
 244       call fitsq(rms,cc(1,1),creff(1,1),ii,przes,obrot,non_conv)
 245       if (non_conv) then
 246         write (iout,*) 'Error: FITSQ non-convergent'
 247         rms=1.0d2
 248       else if (rms.lt.-1.0d-6) then
 249         print *,'Error: rms^2 = ',rms
 250         rms = 1.0d2
 251       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
 252         rmscalc=0.0d0
 253       else
 254         rms = dsqrt(rms)
 255       endif
 256       if (rms.lt.rmsmin) then
 257         rmsmin=rms
 258         ipermmin=iperm
 259       endif
 260 #ifdef DEBUG
 261       write (iout,*) "iperm",iperm," rms",rms
 262 #endif
 263       ENDDO
 264       rmscalc=rmsmin
 265 #ifdef DEBUG
 266       write (iout,*) "ipermmin",ipermmin," rmsmin",rmsmin
 267 #endif
 268       return
 269       end
 270 c-----------------------------------------------------------------------------
 271       double precision function gyrate(jcon)
 272       implicit real*8 (a-h,o-z)
 273       include 'DIMENSIONS'
 274       include 'COMMON.INTERACT'
 275       include 'COMMON.CHAIN'
 276       include 'COMMON.IOUNITS'
 277       double precision cen(3),rg
 278
 279       do j=1,3
 280        cen(j)=0.0d0
 281       enddo
 282
 283       ii=0
 284       do i=nnt,nct
 285         if (itype(i).eq.ntyp1) cycle
 286         ii=ii+1
 287         do j=1,3
 288           cen(j)=cen(j)+c(j,i)
 289         enddo
 290       enddo
 291       do j=1,3
 292         cen(j)=cen(j)/dble(ii)
 293       enddo
 294       rg = 0.0d0
 295       do i = nnt, nct
 296         if (itype(i).eq.ntyp1) cycle
 297         do j=1,3
 298          rg = rg + (c(j,i)-cen(j))**2
 299         enddo
 300       end do
 301       gyrate = dsqrt(rg/dble(ii))
 302       return
 303       end