source/cluster/wham/src-M/wrtclust.f

   1       SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'sizesclu.dat'
   5       parameter (num_in_line=5)
   6       LOGICAL PRINTANG(max_cut)
   7       integer PRINTPDB(max_cut),printmol2(max_cut)
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.HEADER'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.VAR'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.IOUNITS'
  14       include 'COMMON.GEO'
  15       include 'COMMON.FREE'
  16       include 'COMMON.TEMPFAC'
  17       CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
  18       character*80 cfname
  19       character*8 ctemper
  20       DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
  21       external ilen
  22       integer ib
  23
  24       do i=1,64
  25         cfname(i:i)=" "
  26       enddo
  27 c      print *,"calling WRTCLUST",ncon
  28 c      write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
  29       rewind 80
  30       call flush(iout)
  31       temper=1.0d0/(beta_h(ib)*1.987d-3)
  32       if (temper.lt.100.0d0) then
  33         write(ctemper,'(f3.0)') temper
  34         ctemper(3:3)=" "
  35       else if (temper.lt.1000.0) then
  36         write (ctemper,'(f4.0)') temper
  37         ctemper(4:4)=" "
  38       else
  39         write (ctemper,'(f5.0)') temper
  40         ctemper(5:5)=" "
  41       endif
  42
  43       do i=1,ncon*(ncon-1)/2
  44         read (80) diss(i)
  45       enddo
  46       close(80,status='delete')
  47 C
  48 C  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
  49 C
  50       ii1= index(intinname,'/')
  51       ii2=ii1
  52       ii1=ii1+1
  53       do while (ii2.gt.0)
  54         ii1=ii1+ii2
  55         ii2=index(intinname(ii1:),'/')
  56       enddo
  57       ii = ii1+index(intinname(ii1:),'.')-1
  58       if (ii.eq.0) then
  59         ii=ilen(intinname)
  60       else
  61         ii=ii-1
  62       endif
  63       prefixp=intinname(ii1:ii)
  64 cd    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
  65 cd    print *,'ecut=',ecut
  66       WRITE (iout,100) NGR
  67       DO 19 IGR=1,NGR
  68       WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
  69       NRECORD=LICZ(IGR)/num_in_line
  70       IND1=1
  71       DO 63 IRECORD=1,NRECORD
  72       IND2=IND1+num_in_line-1
  73       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
  74      &    totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
  75       IND1=IND2+1
  76    63 CONTINUE
  77       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
  78      &   totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
  79       IND1=1
  80       ICON=list_conf(NCONF(IGR,1))
  81 c      WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
  82 C 12/8/93 Estimation of "diameters" of the subsequent families.
  83       ave_dim=0.0
  84       amax_dim=0.0
  85 c      write (iout,*) "ecut",ecut
  86       emin=totfree(nconf(igr,1))
  87       do i=2,licz(igr)
  88         ii=nconf(igr,i)
  89         if (totfree(ii)-emin .gt. ecut) goto 10
  90         do j=1,i-1
  91           jj=nconf(igr,j)
  92           if (jj.eq.1) exit
  93           if (ii.lt.jj) then
  94             ind=ioffset(ncon,ii,jj)
  95           else
  96             ind=ioffset(ncon,jj,ii)
  97           endif
  98 c          write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
  99 c     &     " ind",ind
 100           call flush(iout)
 101           curr_dist=dabs(diss(ind)+0.0d0)
 102 c          write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
 103 c     &      list_conf(jj),curr_dist
 104           if (curr_dist .gt. amax_dim) amax_dim=curr_dist
 105           ave_dim=ave_dim+curr_dist**2
 106         enddo
 107       enddo
 108    10 if (licz(igr) .gt. 1)
 109      & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
 110       write (iout,'(/A,F8.1,A,F8.1)')
 111      & 'Max. distance in the family:',amax_dim,
 112      & '; average distance in the family:',ave_dim
 113       rmsave(igr)=0.0d0
 114       gdt_ts_ave(igr)=0.0d0
 115       gdt_ha_ave(igr)=0.0d0
 116       tmscore_ave(igr)=0.0d0
 117       qpart=0.0d0
 118       e1=totfree(nconf(igr,1))
 119       do i=1,licz(igr)
 120         icon=nconf(igr,i)
 121         boltz=dexp(-(totfree(icon)-e1))
 122         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
 123         gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
 124         gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
 125         tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
 126         qpart=qpart+boltz
 127       enddo
 128       rmsave(igr)=rmsave(igr)/qpart
 129       gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
 130       gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
 131       tmscore_ave(igr)=tmscore_ave(igr)/qpart
 132       write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
 133      & rmsave(igr)," A, ",
 134      & "TMscore",tmscore_ave(igr),
 135      & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
 136      & gdt_ha_ave(igr)
 137    19 CONTINUE
 138       WRITE (iout,400)
 139       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
 140 c      print *,icut,printang(icut)
 141       IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
 142         emin=totfree_gr(1)
 143 c        print *,'emin',emin,' ngr',ngr
 144         if (lprint_cart) then
 145           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 146      &      //"K"//".x"
 147         else
 148           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 149      &      //"K"//".int"
 150         endif
 151         do igr=1,ngr
 152           icon=nconf(igr,1)
 153           if (totfree_gr(igr)-emin.le.ecut) then
 154             if (lprint_cart) then
 155               call cartout(igr,icon,totfree(icon)/beta_h(ib),
 156      &          totfree_gr(igr)/beta_h(ib),
 157      &          rmstb(icon),cfname)
 158             else
 159 c              print '(a)','calling briefout'
 160               do i=1,2*nres
 161                 do j=1,3
 162                   c(j,i)=allcart(j,i,icon)
 163                 enddo
 164               enddo
 165               call int_from_cart1(.false.)
 166               call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
 167      &          totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
 168      &          jhpb_all(1,icon),cfname)
 169 c              print '(a)','exit briefout'
 170             endif
 171           endif
 172         enddo
 173         close(igeom)
 174       ENDIF
 175       IF (PRINTPDB(ICUT).gt.0) THEN
 176 c Write out a number of conformations from each family in PDB format and
 177 c create InsightII command file for their displaying in different colors
 178         cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 179      &    //"K_"//'ave'//exten
 180         write (iout,*) "cfname",cfname
 181         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 182         write (ipdb,'(a,f8.2)')
 183      &    "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
 184         close (ipdb)
 185         I=1
 186         ICON=NCONF(1,1)
 187         EMIN=totfree_gr(I)
 188         emin1=totfree(icon)
 189         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 190 c          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
 191           write (NUMM,'(bz,i4.4)') i
 192           ncon_lim=min0(licz(i),printpdb(icut))
 193           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 194      &      //"K_"//numm(:ilen(numm))//exten
 195           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 196           write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
 197      &     " AVE RMSD",0pf5.2)')
 198      &     i,totfree_gr(i)/beta_h(ib),rmsave(i)
 199 c Write conformations of the family i to PDB files
 200           ncon_out=1
 201           do while (ncon_out.lt.printpdb(icut) .and.
 202      &     ncon_out.lt.licz(i).and.
 203      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 204             ncon_out=ncon_out+1
 205 c            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 206 c     &        totfree(nconf(i,ncon_out)),emin1,ecut
 207           enddo
 208 c          write (iout,*) "ncon_out",ncon_out
 209           call flush(iout)
 210           do j=1,nres
 211             tempfac(1,j)=5.0d0
 212             tempfac(2,j)=5.0d0
 213           enddo
 214           do j=1,ncon_out
 215             icon=nconf(i,j)
 216             do ii=1,2*nres
 217               do k=1,3
 218                 c(k,ii)=allcart(k,ii,icon)
 219               enddo
 220             enddo
 221             call center
 222             call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 223             write (ipdb,'("TER")')
 224           enddo
 225           close(ipdb)
 226 c Average structures and structures closest to average
 227           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 228      &    //"K_"//'ave'//exten
 229           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
 230      &     position="APPEND")
 231           call ave_coord(i)
 232           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 233           call center
 234           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 235           write (ipdb,'("TER")')
 236           call closest_coord(i)
 237 c          write (iout,*) "Calling rmsnat"
 238           rms_closest(i) = rmsnat(i)
 239           call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
 240      &      tmscore_closest(i),cfname,.true.)
 241 c          write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
 242 c     &      gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
 243 c     &      "tmscore",tmscore_closest(i)
 244           call center
 245           call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
 246           write (ipdb,'("TER")')
 247           close (ipdb)
 248           I=I+1
 249           ICON=NCONF(I,1)
 250           emin1=totfree(icon)
 251         ENDDO
 252       ENDIF
 253       IF (printmol2(icut).gt.0) THEN
 254 c Write out a number of conformations from each family in PDB format and
 255 c create InsightII command file for their displaying in different colors
 256         I=1
 257         ICON=NCONF(1,1)
 258         EMIN=ENERGY(ICON)
 259         emin1=emin
 260         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 261           write (NUMM,'(bz,i4.4)') i
 262           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 263      &      //"K_"//numm(:ilen(numm))//extmol
 264           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 265           ncon_out=1
 266           do while (ncon_out.lt.printmol2(icut) .and.
 267      &     ncon_out.lt.licz(i).and.
 268      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 269             ncon_out=ncon_out+1
 270           enddo
 271           do j=1,ncon_out
 272             icon=nconf(i,j)
 273             do ii=1,2*nres
 274               do k=1,3
 275                 c(k,ii)=allcart(k,ii,icon)
 276               enddo
 277             enddo
 278             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 279           enddo
 280           CLOSE(imol2)
 281           I=I+1
 282           ICON=NCONF(I,1)
 283           emin1=totfree(icon)
 284         ENDDO
 285       ENDIF
 286       call WRITE_STATS(ICUT,NCON,IB)
 287   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 288   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
 289      & ' CONTAINS ',I4,' CONFORMATION(S): ')
 290 c 300 FORMAT ( 8(I4,F6.1))
 291   300 FORMAT (5(I4,1pe12.3))
 292   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 293   500 FORMAT (8(2I4,2X))
 294   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 295       RETURN
 296       END
 297 c------------------------------------------------------------------------------
 298       subroutine ave_coord(igr)
 299       implicit none
 300       include 'DIMENSIONS'
 301       include 'sizesclu.dat'
 302       include 'COMMON.CONTROL'
 303       include 'COMMON.CLUSTER'
 304       include 'COMMON.CHAIN'
 305       include 'COMMON.INTERACT'
 306       include 'COMMON.VAR'
 307       include 'COMMON.TEMPFAC'
 308       include 'COMMON.IOUNITS'
 309       logical non_conv
 310       double precision przes(3),obrot(3,3)
 311       double precision xx(3,maxres2),yy(3,maxres2),csq(3,maxres2)
 312       double precision eref
 313       integer i,ii,j,k,icon,jcon,igr
 314       double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
 315 c      write (iout,*) "AVE_COORD: igr",igr
 316       jcon=nconf(igr,1)
 317       eref=totfree(jcon)
 318       boltz = dexp(-totfree(jcon)+eref)
 319       qpart=boltz
 320       do i=1,2*nres
 321         do j=1,3
 322           c(j,i)=allcart(j,i,jcon)*boltz
 323           cref1(j,i)=allcart(j,i,jcon)
 324           csq(j,i)=allcart(j,i,jcon)**2*boltz
 325         enddo
 326       enddo
 327       DO K=2,LICZ(IGR)
 328       jcon=nconf(igr,k)
 329       if (lside) then
 330         ii=0
 331         do i=nnt,nct
 332           ii=ii+1
 333           do j=1,3
 334             xx(j,ii)=allcart(j,i,jcon)
 335             yy(j,ii)=cref1(j,i)
 336           enddo
 337         enddo
 338         do i=nnt,nct
 339 c          if (itype(i).ne.10) then
 340             ii=ii+1
 341             do j=1,3
 342               xx(j,ii)=allcart(j,i+nres,jcon)
 343               yy(j,ii)=cref1(j,i+nres)
 344             enddo
 345 c          endif
 346         enddo
 347         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 348       else
 349         do i=nnt,nct
 350           do j=1,3
 351             cwork(j,i)=allcart(j,i,jcon)
 352           enddo
 353         enddo
 354         call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot
 355      &       ,non_conv)
 356       endif
 357 c      write (iout,*) "rms",rms
 358 c      do i=1,3
 359 c        write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
 360 c      enddo
 361       if (rms.lt.0.0) then
 362         print *,'error, rms^2 = ',rms,icon,jcon
 363         stop
 364       endif
 365       if (non_conv) print *,non_conv,icon,jcon
 366       boltz=dexp(-totfree(jcon)+eref)
 367       qpart = qpart + boltz
 368       do i=1,2*nres
 369         do j=1,3
 370           xx(j,i)=allcart(j,i,jcon)
 371         enddo
 372       enddo
 373       call matvec(cwork,obrot,xx,2*nres)
 374       do i=1,2*nres
 375 c        write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 376 c     &    (allcart(j,i,jcon),j=1,3)
 377         do j=1,3
 378           cwork(j,i)=cwork(j,i)+przes(j)
 379           c(j,i)=c(j,i)+cwork(j,i)*boltz
 380           csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 381         enddo
 382       enddo
 383       ENDDO ! K
 384       do i=1,2*nres
 385         do j=1,3
 386           c(j,i)=c(j,i)/qpart
 387           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 388         enddo
 389 c        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 390       enddo
 391       do i=nnt,nct
 392         tempfac(1,i)=0.0d0
 393         tempfac(2,i)=0.0d0
 394         do j=1,3
 395           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 396           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 397         enddo
 398         tempfac(1,i)=dsqrt(tempfac(1,i))
 399         tempfac(2,i)=dsqrt(tempfac(2,i))
 400       enddo
 401       return
 402       end
 403 c------------------------------------------------------------------------------
 404       subroutine closest_coord(igr)
 405       implicit none
 406       include 'DIMENSIONS'
 407       include 'sizesclu.dat'
 408       include 'COMMON.IOUNITS'
 409       include 'COMMON.CONTROL'
 410       include 'COMMON.CLUSTER'
 411       include 'COMMON.CHAIN'
 412       include 'COMMON.INTERACT'
 413       include 'COMMON.VAR'
 414       logical non_conv
 415       double precision przes(3),obrot(3,3)
 416       double precision xx(3,maxres2),yy(3,maxres2)
 417       integer i,ii,j,k,icon,jcon,jconmin,igr
 418       double precision rms,rmsmin,cwork(3,maxres2)
 419       rmsmin=1.0d10
 420       jconmin=nconf(igr,1)
 421       DO K=1,LICZ(IGR)
 422       jcon=nconf(igr,k)
 423       if (lside) then
 424         ii=0
 425         do i=nnt,nct
 426           ii=ii+1
 427           do j=1,3
 428             xx(j,ii)=allcart(j,i,jcon)
 429             yy(j,ii)=c(j,i)
 430           enddo
 431         enddo
 432         do i=nnt,nct
 433 c          if (itype(i).ne.10) then
 434             ii=ii+1
 435             do j=1,3
 436               xx(j,ii)=allcart(j,i+nres,jcon)
 437               yy(j,ii)=c(j,i+nres)
 438             enddo
 439 c          endif
 440         enddo
 441         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 442       else
 443         do i=nnt,nct
 444           do j=1,3
 445             cwork(j,i)=allcart(j,i,jcon)
 446           enddo
 447         enddo
 448         call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot
 449      &       ,non_conv)
 450       endif
 451       if (rms.lt.0.0) then
 452         print *,'error, rms^2 = ',rms,icon,jcon
 453         stop
 454       endif
 455 c      write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 456       if (non_conv) print *,non_conv,icon,jcon
 457       if (rms.lt.rmsmin) then
 458         rmsmin=rms
 459         jconmin=jcon
 460       endif
 461       ENDDO ! K
 462 c      write (iout,*) "rmsmin",rmsmin," rms",rms
 463       call flush(iout)
 464       do i=1,2*nres
 465         do j=1,3
 466           c(j,i)=allcart(j,i,jconmin)
 467         enddo
 468       enddo
 469       return
 470       end
 471 c------------------------------------------------------------------------------
 472       subroutine center
 473       implicit none
 474       include 'DIMENSIONS'
 475       include 'sizesclu.dat'
 476       include 'COMMON.IOUNITS'
 477       include 'COMMON.CONTROL'
 478       include 'COMMON.CLUSTER'
 479       include 'COMMON.CHAIN'
 480       include 'COMMON.INTERACT'
 481       include 'COMMON.VAR'
 482       double precision przes(3)
 483       integer i,ii,j,k,icon,jcon,jconmin,igr
 484       przes=0.0d0
 485       do j=1,3
 486         do i=1,nres
 487           przes(j)=przes(j)+c(j,i)
 488         enddo
 489       enddo
 490       do j=1,3
 491         przes(j)=przes(j)/nres
 492       enddo
 493       do i=1,2*nres
 494         do j=1,3
 495           c(j,i)=c(j,i)-przes(j)
 496         enddo
 497       enddo
 498       return
 499       end