source/cluster/unres/src/wrtclust.f

   1       SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2)
   2       include 'DIMENSIONS'
   3       include 'sizesclu.dat'
   4       parameter (num_in_line=5)
   5       LOGICAL PRINTANG(max_cut),linsight
   6       integer PRINTPDB(max_cut),printmol2(max_cut)
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.HEADER'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CLUSTER'
  12       include 'COMMON.IOUNITS'
  13       include 'COMMON.GEO'
  14       double precision totfree_gr(maxconf)
  15       CHARACTER*64 prefixp,CFNAME,CFNAME1,CFF,NUMM,MUMM,EXTEN,UCASE,
  16      &  extmol
  17       DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
  18       external ilen
  19       logical insight_cmd_out
  20
  21 c      print *,"calling WRTCLUST",ncon
  22 c      ICANT(I,J)=((NCON+NCON-J)*(J-1))/2+I-J
  23       insight_cmd_out = .false.
  24
  25 C
  26 C  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
  27 C
  28       ii1= index(intinname,'/')
  29       ii2=ii1
  30       ii1=ii1+1
  31       do while (ii2.gt.0)
  32         ii1=ii1+ii2
  33         ii2=index(intinname(ii1:),'/')
  34       enddo
  35       ii = ii1+index(intinname(ii1:),'.')-1
  36       if (ii.eq.0) then
  37         ii=ilen(intinname)
  38       else
  39         ii=ii-1
  40       endif
  41       prefixp=intinname(ii1:ii)
  42 cd    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
  43 cd    print *,'ecut=',ecut
  44       WRITE (iout,100) NGR
  45       DO 19 IGR=1,NGR
  46       WRITE (iout,200) IGR,LICZ(IGR)
  47       NRECORD=LICZ(IGR)/num_in_line
  48       IND1=1
  49       DO 63 IRECORD=1,NRECORD
  50       IND2=IND1+num_in_line-1
  51       WRITE (iout,300) (NCONF(IGR,ICO),ENERGY(NCONF(IGR,ICO)),
  52      1              ICO=IND1,IND2)
  53       IND1=IND2+1
  54    63 CONTINUE
  55       WRITE (iout,300) (NCONF(IGR,ICO),ENERGY(NCONF(IGR,ICO)),
  56      1              ICO=IND1,LICZ(IGR))
  57       IND1=1
  58       ICON=NCONF(IGR,1)
  59       WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
  60 C 12/8/93 Estimation of "diameters" of the subsequent families.
  61       emin=totfree(nconf(igr,1))
  62       if (efree) then
  63         totfree_gr(igr)=1.0d0
  64         do i=2,licz(igr)
  65           ii=nconf(igr,i)
  66           totfree_gr(igr)=totfree_gr(igr)
  67      &     +dexp(-betaT*(totfree(ii)-emin))
  68         enddo
  69         write (iout,*) "igr",igr," totfree",emin,
  70      &    " totfree_gr",totfree_gr(igr)
  71         totfree_gr(igr)=emin-dlog(totfree_gr(igr))/betaT
  72         write (iout,*) "efree",totfree_gr(igr)
  73       endif
  74       ave_dim=0.0
  75       amax_dim=0.0
  76       do i=2,licz(igr)
  77         ii=nconf(igr,i)
  78         if (totfree(ii)-emin .gt. ecut) goto 10
  79         do j=1,i-1
  80           jj=nconf(igr,j)
  81           if (ii.lt.jj) then
  82             ind=ioffset(ncon,ii,jj)
  83           else
  84             ind=ioffset(ncon,jj,ii)
  85           endif
  86           curr_dist=dabs(attalums(ind))
  87 cd          print '(3i4,f12.4)',ind,ii,jj,curr_dist
  88           if (curr_dist .gt. amax_dim) amax_dim=curr_dist
  89           ave_dim=ave_dim+curr_dist**2
  90         enddo
  91       enddo
  92    10 if (licz(igr) .gt. 1)
  93      & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
  94       write (iout,'(/A,F8.1,A,F8.1)')
  95      & 'Max. distance in the family:',amax_dim,
  96      & '; average distance in the family:',ave_dim
  97       if (refstr .or. pdbref) write (iout,'(a,i5,f8.3)')
  98      & "RMSD of the lowest-energy conformation #",nconf(igr,1),
  99      &  rmsnat(nconf(igr,1))
 100    19 CONTINUE
 101       WRITE (iout,400)
 102       WRITE (iout,500) (I,IASS(I),I=1,NCON)
 103 c      print *,icut,printang(icut)
 104       IF (PRINTANG(ICUT)) then
 105         emin=totfree(nconf(1,1))
 106 c        print *,'emin',emin,' ngr',ngr
 107         do igr=1,ngr
 108           icon=nconf(igr,1)
 109           if (totfree(icon)-emin.le.ecut) then
 110             do i=1,nres
 111               phi(i)=phiall(i,icon)
 112               theta(i)=thetall(i,icon)
 113               alph(i)=alphall(i,icon)
 114               omeg(i)=omall(i,icon)
 115             enddo
 116             if (lprint_cart) then
 117               call cartout(igr,icon,energy(icon),totfree_gr(igr),
 118      &          rmstab(icon),intname)
 119             else
 120 c              print '(a)','calling briefout'
 121               call briefout(igr,iscore(icon),energy(icon),
 122      &          totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
 123      &          jhpb_all(1,icon),intname)
 124 c              print '(a)','exit briefout'
 125             endif
 126           endif
 127         enddo
 128       ENDIF
 129       IF (PRINTPDB(ICUT).gt.0) THEN
 130 c Write out a number of conformations from each family in PDB format and
 131 c create InsightII command file for their displaying in different colors
 132         I=1
 133         ICON=NCONF(1,1)
 134         EMIN=ENERGY(ICON)
 135         DO WHILE(I.LE.NGR .AND. ENERGY(ICON)-EMIN.LE.ECUT)
 136 c          CALL NUMSTR(I,NUMM)
 137           write (NUMM,'(bz,i4.4)') i
 138           ncon_lim=min0(licz(i),printpdb(icut))
 139           linsight=
 140      &       printpdb(icut).gt.0 .and. ncon_lim.gt.1
 141           if (linsight) then
 142             icon=nconf(i,1)
 143             emini=energy(icon)
 144             k=1
 145             do while (k.le.ncon_lim .and.
 146      &                        energy(nconf(i,k))-emini.le.ecut)
 147               k=k+1
 148             enddo
 149             ncon_out=min0(k-1,ncon_lim)
 150             linsight=ncon_out.gt.1
 151           endif
 152
 153           if (linsight) then
 154 c
 155 c A "bunch of structures" of the family that lie within ECUT above the
 156 c lowest-energy conformation in the family will be outputed along with the
 157 c InsightII command file --- AL 1/1/95.
 158 c
 159             if (insight_cmd_out) then
 160             open (22,file=
 161      & 'insight_'//prefixp(:ilen(prefixp))//numm(:ilen(numm))//'.cmd',
 162      &   status='unknown')
 163 c
 164 c Write InsightII command file
 165 c
 166             cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_001'
 167             cff=ucase(cfname)
 168             write (22,'(5a)') 'get molecule pdb user ',
 169      &            cfname(:ilen(cfname)),'.pdb ',
 170      &            cff(:ilen(cff)),
 171      &            ' -heteroatom -reference_object'
 172             cfname=ucase(cfname)
 173             write (22,'(3a,i3,2(1h,,i3))') 'color molecule atoms ',
 174      &            cfname(:ilen(cfname)),' specified specification ',
 175      &            0,0,255
 176             write (22,'(2a)') 'display molecule only atoms heavy ',
 177      &            cfname(:ilen(cfname))
 178             cfname1=cfname
 179             deltae_max=energy(nconf(i,ncon_out))-emini
 180             do j=2,ncon_out
 181               deltae=energy(nconf(i,j))-emini
 182               icolor=aint(255*deltae/deltae_max)
 183               call numstr(j,mumm)
 184               cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_'
 185      &            //mumm(:ilen(mumm))
 186               cff=ucase(cfname)
 187               write (22,'(5a)') 'get molecule pdb user ',
 188      &            cfname(:ilen(cfname)),'.pdb ',
 189      &            cff(:ilen(cff)),' -heteroatom'
 190               cfname=ucase(cfname)
 191               do k=1,nct-nnt+1
 192                if (k.lt.10) then
 193                 write (22,'(3a,i1,3a,i1)')
 194      &        'superimpose -end_definition backbone "label mode" ',
 195      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 196      &         ':',k
 197                elseif (k.lt.100) then
 198                 write (22,'(3a,i2,3a,i2)')
 199      &        'superimpose -end_definition backbone "label mode" ',
 200      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 201      &         ':',k
 202                else
 203                 write (22,'(3a,i3,3a,i3)')
 204      &        'superimpose -end_definition backbone "label mode" ',
 205      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 206      &         ':',k
 207                endif
 208
 209               enddo
 210               write (22,'(3a,i3,2(1h,,i3))') 'color molecule atoms ',
 211      &            cfname(:ilen(cfname)),' specified specification ',
 212      &            icolor,icolor,255-icolor
 213               write (22,'(2a)') 'display molecule only atoms heavy ',
 214      &            cfname(:ilen(cfname))
 215               close(22)
 216             enddo
 217             endif
 218
 219 c Write conformations of the family i to PDB files
 220             do j=1,ncon_out
 221               call numstr(j,mumm)
 222               icon=nconf(i,j)
 223               cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_'
 224      &             //mumm(:ilen(mumm))//exten
 225               do ii=1,nres
 226                 phi(ii)=phiall(ii,icon)
 227                 theta(ii)=thetall(ii,icon)
 228                 alph(ii)=alphall(ii,icon)
 229                 omeg(ii)=omall(ii,icon)
 230               enddo
 231               call chainbuild
 232               if (refstr.or.pdbref) rmsd=rmsnat(icon)
 233               open(ipdb,file=cfname,status='unknown',form='formatted')
 234               call pdbout(energy(icon),rmsd,titel)
 235               close(ipdb)
 236             enddo
 237           else
 238 c Produce only a single PDB file for the leading member of the family
 239             write (iout,*) 'Writing pdb file: icon=',icon
 240             do ii=1,nres
 241               phi(ii)=phiall(ii,icon)
 242               theta(ii)=thetall(ii,icon)
 243               alph(ii)=alphall(ii,icon)
 244               omeg(ii)=omall(ii,icon)
 245             enddo
 246             call chainbuild
 247             cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//exten
 248             OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 249 c           print *,'Calling pdbout'
 250             if (refstr.or.pdbref) rmsd=rmsnat(icon)
 251             CALL PDBOUT(energy(icon),rmsd,titel)
 252             CLOSE(ipdb)
 253           endif
 254           I=I+1
 255           ICON=NCONF(I,1)
 256         ENDDO
 257       ENDIF
 258       IF (printmol2(icut).gt.0) THEN
 259 c Write out a number of conformations from each family in PDB format and
 260 c create InsightII command file for their displaying in different colors
 261         I=1
 262         ICON=NCONF(1,1)
 263         EMIN=ENERGY(ICON)
 264         DO WHILE(I.LE.NGR .AND. ENERGY(ICON)-EMIN.LE.ECUT)
 265 c          CALL NUMSTR(I,NUMM)
 266           write (NUMM,'(bz,i4.4)') i
 267           ncon_lim=min0(licz(i),printpdb(icut))
 268           write (iout,*) 'Writing mol2 file: icon=',icon
 269           do ii=1,nres
 270             phi(ii)=phiall(ii,icon)
 271             theta(ii)=thetall(ii,icon)
 272             alph(ii)=alphall(ii,icon)
 273             omeg(ii)=omall(ii,icon)
 274           enddo
 275           call chainbuild
 276           cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//extmol
 277           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 278 c         print *,'Calling pdbout'
 279           CALL MOL2OUT(energy(icon),'STRUCTURE'//numm)
 280           CLOSE(imol2)
 281           I=I+1
 282           ICON=NCONF(I,1)
 283         ENDDO
 284       ENDIF
 285   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 286   200 FORMAT (/'FAMILY ',I4,' CONTAINS ',I4,' CONFORMATION(S): ')
 287 c 300 FORMAT ( 8(I4,F6.1))
 288   300 FORMAT (5(I4,1pe12.4))
 289   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 290   500 FORMAT (8(2I4,2X))
 291   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 292       RETURN
 293       END