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