update
[unres.git] / source / cluster / wham / src-M / wrtclust.f
index f2f3eb7..5460eb6 100644 (file)
@@ -14,7 +14,6 @@
       include 'COMMON.GEO'
       include 'COMMON.FREE'
       include 'COMMON.TEMPFAC'
-      double precision rmsave(maxgr)
       CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
       character*80 cfname
       character*8 ctemper
@@ -84,6 +83,7 @@ C 12/8/93 Estimation of "diameters" of the subsequent families.
       ave_dim=0.0
       amax_dim=0.0
 c      write (iout,*) "ecut",ecut
+      emin=totfree(nconf(igr,1))
       do i=2,licz(igr)
         ii=nconf(igr,i)
         if (totfree(ii)-emin .gt. ecut) goto 10
@@ -111,15 +111,29 @@ c     &      list_conf(jj),curr_dist
      & 'Max. distance in the family:',amax_dim,
      & '; average distance in the family:',ave_dim 
       rmsave(igr)=0.0d0
+      gdt_ts_ave(igr)=0.0d0
+      gdt_ha_ave(igr)=0.0d0
+      tmscore_ave(igr)=0.0d0
       qpart=0.0d0
+      e1=totfree(nconf(igr,1))
       do i=1,licz(igr)
         icon=nconf(igr,i)
-        boltz=dexp(-totfree(icon))
+        boltz=dexp(-(totfree(icon)-e1))
         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
+        gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
+        gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
+        tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
         qpart=qpart+boltz
       enddo
       rmsave(igr)=rmsave(igr)/qpart
-      write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
+      gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
+      gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
+      tmscore_ave(igr)=tmscore_ave(igr)/qpart
+      write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
+     & rmsave(igr)," A, ",
+     & "TMscore",tmscore_ave(igr),
+     & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
+     & gdt_ha_ave(igr)
    19 CONTINUE
       WRITE (iout,400)
       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
@@ -191,7 +205,7 @@ c Write conformations of the family i to PDB files
 c            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 c     &        totfree(nconf(i,ncon_out)),emin1,ecut
           enddo
-          write (iout,*) "ncon_out",ncon_out
+c          write (iout,*) "ncon_out",ncon_out
           call flush(iout)
           do j=1,nres
             tempfac(1,j)=5.0d0
@@ -204,6 +218,7 @@ c     &        totfree(nconf(i,ncon_out)),emin1,ecut
                 c(k,ii)=allcart(k,ii,icon)
               enddo
             enddo
+            call center
             call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
             write (ipdb,'("TER")')
           enddo
@@ -215,10 +230,19 @@ c Average structures and structures closest to average
      &     position="APPEND")
           call ave_coord(i)
           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
+          call center
           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
           write (ipdb,'("TER")')
           call closest_coord(i)
-          call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
+c          write (iout,*) "Calling rmsnat"
+          rms_closest(i) = rmsnat(i)
+          call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
+     &      tmscore_closest(i),cfname,.true.)
+c          write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
+c     &      gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
+c     &      "tmscore",tmscore_closest(i)
+          call center
+          call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
           write (ipdb,'("TER")')
           close (ipdb)
           I=I+1
@@ -259,6 +283,7 @@ c create InsightII command file for their displaying in different colors
           emin1=totfree(icon)
         ENDDO
       ENDIF 
+      call WRITE_STATS(ICUT,NCON,IB)
   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
      & ' CONTAINS ',I4,' CONFORMATION(S): ')
@@ -443,3 +468,32 @@ c      write (iout,*) "rmsmin",rmsmin," rms",rms
       enddo
       return
       end
+c------------------------------------------------------------------------------
+      subroutine center
+      implicit none
+      include 'DIMENSIONS'
+      include 'sizesclu.dat'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CLUSTER'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      double precision przes(3)
+      integer i,ii,j,k,icon,jcon,jconmin,igr
+      przes=0.0d0
+      do j=1,3
+        do i=1,nres
+          przes(j)=przes(j)+c(j,i)
+        enddo
+      enddo
+      do j=1,3
+        przes(j)=przes(j)/nres
+      enddo
+      do i=1,2*nres
+        do j=1,3
+          c(j,i)=c(j,i)-przes(j)  
+        enddo
+      enddo
+      return
+      end