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       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       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       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       qpart=0.0d0
 115       do i=1,licz(igr)
 116         icon=nconf(igr,i)
 117         boltz=dexp(-totfree(icon))
 118         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
 119         qpart=qpart+boltz
 120       enddo
 121       rmsave(igr)=rmsave(igr)/qpart
 122       write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
 123    19 CONTINUE
 124       WRITE (iout,400)
 125       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
 126 c      print *,icut,printang(icut)
 127       IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
 128         emin=totfree_gr(1)
 129 c        print *,'emin',emin,' ngr',ngr
 130         if (lprint_cart) then
 131           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 132      &      //"K"//".x"
 133         else
 134           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 135      &      //"K"//".int"
 136         endif
 137         do igr=1,ngr
 138           icon=nconf(igr,1)
 139           if (totfree_gr(igr)-emin.le.ecut) then
 140             if (lprint_cart) then
 141               call cartout(igr,icon,totfree(icon)/beta_h(ib),
 142      &          totfree_gr(igr)/beta_h(ib),
 143      &          rmstb(icon),cfname)
 144             else
 145 c              print '(a)','calling briefout'
 146               do i=1,2*nres
 147                 do j=1,3
 148                   c(j,i)=allcart(j,i,icon)
 149                 enddo
 150               enddo
 151               call int_from_cart1(.false.)
 152               call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
 153      &          totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
 154      &          jhpb_all(1,icon),cfname)
 155 c              print '(a)','exit briefout'
 156             endif
 157           endif
 158         enddo
 159         close(igeom)
 160       ENDIF
 161       IF (PRINTPDB(ICUT).gt.0) THEN
 162 c Write out a number of conformations from each family in PDB format and
 163 c create InsightII command file for their displaying in different colors
 164         cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 165      &    //"K_"//'ave'//exten
 166         write (iout,*) "cfname",cfname
 167         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 168         write (ipdb,'(a,f8.2)')
 169      &    "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
 170         close (ipdb)
 171         I=1
 172         ICON=NCONF(1,1)
 173         EMIN=totfree_gr(I)
 174         emin1=totfree(icon)
 175         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 176 c          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
 177           write (NUMM,'(bz,i4.4)') i
 178           ncon_lim=min0(licz(i),printpdb(icut))
 179           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 180      &      //"K_"//numm(:ilen(numm))//exten
 181           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 182           write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
 183      &     " AVE RMSD",0pf5.2)')
 184      &     i,totfree_gr(i)/beta_h(ib),rmsave(i)
 185 c Write conformations of the family i to PDB files
 186           ncon_out=1
 187           do while (ncon_out.lt.printpdb(icut) .and.
 188      &     ncon_out.lt.licz(i).and.
 189      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 190             ncon_out=ncon_out+1
 191 c            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 192 c     &        totfree(nconf(i,ncon_out)),emin1,ecut
 193           enddo
 194           write (iout,*) "ncon_out",ncon_out
 195           call flush(iout)
 196           do j=1,nres
 197             tempfac(1,j)=5.0d0
 198             tempfac(2,j)=5.0d0
 199           enddo
 200           do j=1,ncon_out
 201             icon=nconf(i,j)
 202             do ii=1,2*nres
 203               do k=1,3
 204                 c(k,ii)=allcart(k,ii,icon)
 205               enddo
 206             enddo
 207             call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 208             write (ipdb,'("TER")')
 209           enddo
 210           close(ipdb)
 211 c Average structures and structures closest to average
 212           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 213      &    //"K_"//'ave'//exten
 214           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
 215      &     position="APPEND")
 216           call ave_coord(i)
 217           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 218           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 219           write (ipdb,'("TER")')
 220           call closest_coord(i)
 221           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 222           write (ipdb,'("TER")')
 223           close (ipdb)
 224           I=I+1
 225           ICON=NCONF(I,1)
 226           emin1=totfree(icon)
 227         ENDDO
 228       ENDIF
 229       IF (printmol2(icut).gt.0) THEN
 230 c Write out a number of conformations from each family in PDB format and
 231 c create InsightII command file for their displaying in different colors
 232         I=1
 233         ICON=NCONF(1,1)
 234         EMIN=ENERGY(ICON)
 235         emin1=emin
 236         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 237           write (NUMM,'(bz,i4.4)') i
 238           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 239      &      //"K_"//numm(:ilen(numm))//extmol
 240           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 241           ncon_out=1
 242           do while (ncon_out.lt.printmol2(icut) .and.
 243      &     ncon_out.lt.licz(i).and.
 244      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 245             ncon_out=ncon_out+1
 246           enddo
 247           do j=1,ncon_out
 248             icon=nconf(i,j)
 249             do ii=1,2*nres
 250               do k=1,3
 251                 c(k,ii)=allcart(k,ii,icon)
 252               enddo
 253             enddo
 254             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 255           enddo
 256           CLOSE(imol2)
 257           I=I+1
 258           ICON=NCONF(I,1)
 259           emin1=totfree(icon)
 260         ENDDO
 261       ENDIF
 262   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 263   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
 264      & ' CONTAINS ',I4,' CONFORMATION(S): ')
 265 c 300 FORMAT ( 8(I4,F6.1))
 266   300 FORMAT (5(I4,1pe12.3))
 267   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 268   500 FORMAT (8(2I4,2X))
 269   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 270       RETURN
 271       END
 272 c------------------------------------------------------------------------------
 273       subroutine ave_coord(igr)
 274       implicit none
 275       include 'DIMENSIONS'
 276       include 'sizesclu.dat'
 277       include 'COMMON.CONTROL'
 278       include 'COMMON.CLUSTER'
 279       include 'COMMON.CHAIN'
 280       include 'COMMON.INTERACT'
 281       include 'COMMON.VAR'
 282       include 'COMMON.TEMPFAC'
 283       include 'COMMON.IOUNITS'
 284       logical non_conv
 285       double precision przes(3),obrot(3,3)
 286       double precision xx(3,maxres2),yy(3,maxres2),csq(3,maxres2)
 287       double precision eref
 288       integer i,ii,j,k,icon,jcon,igr
 289       double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
 290 c      write (iout,*) "AVE_COORD: igr",igr
 291       jcon=nconf(igr,1)
 292       eref=totfree(jcon)
 293       boltz = dexp(-totfree(jcon)+eref)
 294       qpart=boltz
 295       do i=1,2*nres
 296         do j=1,3
 297           c(j,i)=allcart(j,i,jcon)*boltz
 298           cref1(j,i)=allcart(j,i,jcon)
 299           csq(j,i)=allcart(j,i,jcon)**2*boltz
 300         enddo
 301       enddo
 302       DO K=2,LICZ(IGR)
 303       jcon=nconf(igr,k)
 304       if (lside) then
 305         ii=0
 306         do i=nnt,nct
 307           ii=ii+1
 308           do j=1,3
 309             xx(j,ii)=allcart(j,i,jcon)
 310             yy(j,ii)=cref1(j,i)
 311           enddo
 312         enddo
 313         do i=nnt,nct
 314 c          if (itype(i).ne.10) then
 315             ii=ii+1
 316             do j=1,3
 317               xx(j,ii)=allcart(j,i+nres,jcon)
 318               yy(j,ii)=cref1(j,i+nres)
 319             enddo
 320 c          endif
 321         enddo
 322         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 323       else
 324         do i=nnt,nct
 325           do j=1,3
 326             cwork(j,i)=allcart(j,i,jcon)
 327           enddo
 328         enddo
 329         call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot
 330      &       ,non_conv)
 331       endif
 332 c      write (iout,*) "rms",rms
 333 c      do i=1,3
 334 c        write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
 335 c      enddo
 336       if (rms.lt.0.0) then
 337         print *,'error, rms^2 = ',rms,icon,jcon
 338         stop
 339       endif
 340       if (non_conv) print *,non_conv,icon,jcon
 341       boltz=dexp(-totfree(jcon)+eref)
 342       qpart = qpart + boltz
 343       do i=1,2*nres
 344         do j=1,3
 345           xx(j,i)=allcart(j,i,jcon)
 346         enddo
 347       enddo
 348       call matvec(cwork,obrot,xx,2*nres)
 349       do i=1,2*nres
 350 c        write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 351 c     &    (allcart(j,i,jcon),j=1,3)
 352         do j=1,3
 353           cwork(j,i)=cwork(j,i)+przes(j)
 354           c(j,i)=c(j,i)+cwork(j,i)*boltz
 355           csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 356         enddo
 357       enddo
 358       ENDDO ! K
 359       do i=1,2*nres
 360         do j=1,3
 361           c(j,i)=c(j,i)/qpart
 362           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 363         enddo
 364 c        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 365       enddo
 366       do i=nnt,nct
 367         tempfac(1,i)=0.0d0
 368         tempfac(2,i)=0.0d0
 369         do j=1,3
 370           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 371           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 372         enddo
 373         tempfac(1,i)=dsqrt(tempfac(1,i))
 374         tempfac(2,i)=dsqrt(tempfac(2,i))
 375       enddo
 376       return
 377       end
 378 c------------------------------------------------------------------------------
 379       subroutine closest_coord(igr)
 380       implicit none
 381       include 'DIMENSIONS'
 382       include 'sizesclu.dat'
 383       include 'COMMON.IOUNITS'
 384       include 'COMMON.CONTROL'
 385       include 'COMMON.CLUSTER'
 386       include 'COMMON.CHAIN'
 387       include 'COMMON.INTERACT'
 388       include 'COMMON.VAR'
 389       logical non_conv
 390       double precision przes(3),obrot(3,3)
 391       double precision xx(3,maxres2),yy(3,maxres2)
 392       integer i,ii,j,k,icon,jcon,jconmin,igr
 393       double precision rms,rmsmin,cwork(3,maxres2)
 394       rmsmin=1.0d10
 395       jconmin=nconf(igr,1)
 396       DO K=1,LICZ(IGR)
 397       jcon=nconf(igr,k)
 398       if (lside) then
 399         ii=0
 400         do i=nnt,nct
 401           ii=ii+1
 402           do j=1,3
 403             xx(j,ii)=allcart(j,i,jcon)
 404             yy(j,ii)=c(j,i)
 405           enddo
 406         enddo
 407         do i=nnt,nct
 408 c          if (itype(i).ne.10) then
 409             ii=ii+1
 410             do j=1,3
 411               xx(j,ii)=allcart(j,i+nres,jcon)
 412               yy(j,ii)=c(j,i+nres)
 413             enddo
 414 c          endif
 415         enddo
 416         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 417       else
 418         do i=nnt,nct
 419           do j=1,3
 420             cwork(j,i)=allcart(j,i,jcon)
 421           enddo
 422         enddo
 423         call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot
 424      &       ,non_conv)
 425       endif
 426       if (rms.lt.0.0) then
 427         print *,'error, rms^2 = ',rms,icon,jcon
 428         stop
 429       endif
 430 c      write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 431       if (non_conv) print *,non_conv,icon,jcon
 432       if (rms.lt.rmsmin) then
 433         rmsmin=rms
 434         jconmin=jcon
 435       endif
 436       ENDDO ! K
 437 c      write (iout,*) "rmsmin",rmsmin," rms",rms
 438       call flush(iout)
 439       do i=1,2*nres
 440         do j=1,3
 441           c(j,i)=allcart(j,i,jconmin)
 442         enddo
 443       enddo
 444       return
 445       end