Merge branch 'devel' into AFM
[unres.git] / source / wham / src-NEWSC-NEWCORR / rmscalc.f
diff --git a/source/wham/src-NEWSC-NEWCORR/rmscalc.f b/source/wham/src-NEWSC-NEWCORR/rmscalc.f
new file mode 100644 (file)
index 0000000..70d9425
--- /dev/null
@@ -0,0 +1,156 @@
+      double precision function rmscalc(ishif,i,j,jcon,lprn)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.COMPAR'
+      include 'COMMON.CHAIN' 
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      double precision przes(3),obrot(3,3)
+      double precision creff(3,maxres2),cc(3,maxres2)
+      logical iadded(maxres)
+      integer inumber(2,maxres)
+      common /ccc/ creff,cc,iadded,inumber
+      logical lprn
+      logical non_conv
+      integer ishif,i,j
+      if (lprn) then
+      write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
+      write (iout,*) "npiece",npiece(j,i)
+      endif
+      ii=0
+      do l=1,nres
+        iadded(l)=.false.
+      enddo
+      do k=1,npiece(j,i)
+        if (i.eq.1) then
+          if (lprn)
+     &    write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",
+     &     ifrag(1,k,j),ifrag(2,k,j)
+          call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,ii)
+c          write (iout,*) "ii=",ii
+        else
+          kk = ipiece(k,j,i)
+c          write (iout,*) "kk",kk," npiece",npiece(kk,1)
+          do l=1,npiece(kk,1)
+            if (lprn)
+     &      write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,
+     &       " l=",l," adding fragment",
+     &       ifrag(1,l,kk),ifrag(2,l,kk)
+            call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,ii)
+          enddo 
+        endif
+      enddo
+      if (lprn) then
+      do k=1,ii
+        write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),
+     &         inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
+      enddo
+      endif
+      call fitsq(rms,cc(1,1),creff(1,1),ii,przes,obrot,non_conv)
+      if (non_conv) then
+        print *,'Error: FITSQ non-convergent, jcon',jcon
+        rmscalc=1.0d2
+      else if (rms.lt.-1.0d-6) then 
+        print *,'Error: rms^2 = ',rms,jcon
+        rmscalc = 1.0d2
+      else if (rms.ge.1.0d-6 .and. rms.lt.0) then
+        rmscalc=0.0d0
+      else 
+        rmscalc = dsqrt(rms)
+      endif
+      return
+      end
+c-------------------------------------------------------------------------
+      subroutine cprep(if1,if2,ishif,ii)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.COMPAR'
+      include 'COMMON.CHAIN' 
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      double precision przes(3),obrot(3,3)
+      double precision creff(3,maxres2),cc(3,maxres2)
+      logical iadded(maxres)
+      integer inumber(2,maxres)
+      common /ccc/ creff,cc,iadded,inumber
+c      write (iout,*) "Calling cprep"
+      do l=if1,if2
+c        write (iout,*) "l",l," iadded",iadded(l)
+        if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) 
+     &  then
+          ii=ii+1
+          iadded(l)=.true.
+          inumber(1,ii)=l
+          inumber(2,ii)=l+ishif
+          do m=1,3
+            creff(m,ii)=cref(m,l)
+            cc(m,ii)=c(m,l+ishif)
+          enddo
+        endif
+      enddo
+      return
+      end
+c-------------------------------------------------------------------------
+      double precision function rmsnat(jcon)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.COMPAR'
+      include 'COMMON.CHAIN' 
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      double precision przes(3),obrot(3,3)
+      logical non_conv
+      integer ishif,i,j
+      call fitsq(rms,c(1,nstart_sup),cref(1,nstart_sup),nsup,
+     &  przes,obrot,non_conv)
+      if (non_conv) then
+        print *,'Error: FITSQ non-convergent, jcon',jcon
+        rmsnat=1.0d2
+      else if (rms.lt.-1.0d-6) then 
+        print *,'Error: rms^2 = ',rms,jcon
+        rmsnat = 1.0d2
+      else if (rms.ge.1.0d-6 .and. rms.lt.0) then
+        rmsnat=0.0d0
+      else 
+        rmsnat = dsqrt(rms)
+      endif
+      return
+      end
+c-----------------------------------------------------------------------------
+      double precision function gyrate(jcon)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.INTERACT'
+      include 'COMMON.CHAIN'
+      double precision cen(3),rg
+
+      do j=1,3
+       cen(j)=0.0d0
+      enddo
+
+      do i=nnt,nct
+          do j=1,3
+            cen(j)=cen(j)+c(j,i)
+          enddo
+      enddo
+      do j=1,3
+            cen(j)=cen(j)/dble(nct-nnt+1)
+      enddo
+      rg = 0.0d0
+      do i = nnt, nct
+        do j=1,3
+         rg = rg + (c(j,i)-cen(j))**2
+        enddo
+      end do
+      gyrate = dsqrt(rg/dble(nct-nnt+1))
+      return
+      end