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

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