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(intname,'/')
  29       ii2=ii1
  30       ii1=ii1+1
  31       do while (ii2.gt.0)
  32         ii1=ii1+ii2
  33         ii2=index(intname(ii1:),'/')
  34       enddo
  35       ii = ii1+index(intname(ii1:),'.')-1
  36       if (ii.eq.0) then
  37         ii=ilen(intname)
  38       else
  39         ii=ii-1
  40       endif
  41       prefixp=intname(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) then
  98         write (iout,'(a,i5,2f8.3)')
  99      & "RMSD of the lowest-energy conformation #",nconf(igr,1),
 100      &  rmsnat(nconf(igr,1)),rmstab(nconf(igr,1))
 101         rmsave=0.0d0
 102         do i=1,licz(igr)
 103           rmsave=rmsave+rmsnat(nconf(igr,i))
 104         enddo
 105         rmsave=rmsave/licz(igr)
 106         write (iout,'(a,f8.3)') "Average RMSD in the family",rmsave
 107       endif
 108    19 CONTINUE
 109       WRITE (iout,400)
 110       WRITE (iout,500) (I,IASS(I),I=1,NCON)
 111 c      print *,icut,printang(icut)
 112       IF (PRINTANG(ICUT)) then
 113         emin=totfree(nconf(1,1))
 114 c        print *,'emin',emin,' ngr',ngr
 115         do igr=1,ngr
 116           icon=nconf(igr,1)
 117           if (totfree(icon)-emin.le.ecut) then
 118             do i=1,nres
 119               phi(i)=phiall(i,icon)
 120               theta(i)=thetall(i,icon)
 121               alph(i)=alphall(i,icon)
 122               omeg(i)=omall(i,icon)
 123             enddo
 124             if (lprint_cart) then
 125               call cartout(igr,icon,energy(icon),totfree_gr(igr),
 126      &          rmstab(icon),intname)
 127             else
 128 c              print '(a)','calling briefout'
 129               call briefout(igr,iscore(icon),energy(icon),
 130      &          totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
 131      &          jhpb_all(1,icon),intname)
 132 c              print '(a)','exit briefout'
 133             endif
 134           endif
 135         enddo
 136       ENDIF
 137       IF (PRINTPDB(ICUT).gt.0) THEN
 138 c Write out a number of conformations from each family in PDB format and
 139 c create InsightII command file for their displaying in different colors
 140         I=1
 141         ICON=NCONF(1,1)
 142         EMIN=ENERGY(ICON)
 143         DO WHILE(I.LE.NGR .AND. ENERGY(ICON)-EMIN.LE.ECUT)
 144 c          CALL NUMSTR(I,NUMM)
 145           write (NUMM,'(bz,i4.4)') i
 146           ncon_lim=min0(licz(i),printpdb(icut))
 147           linsight=
 148      &       printpdb(icut).gt.0 .and. ncon_lim.gt.1
 149           if (linsight) then
 150             icon=nconf(i,1)
 151             emini=energy(icon)
 152             k=1
 153             do while (k.le.ncon_lim .and.
 154      &                        energy(nconf(i,k))-emini.le.ecut)
 155               k=k+1
 156             enddo
 157             ncon_out=min0(k-1,ncon_lim)
 158             linsight=ncon_out.gt.1
 159           endif
 160
 161           if (linsight) then
 162 c
 163 c A "bunch of structures" of the family that lie within ECUT above the
 164 c lowest-energy conformation in the family will be outputed along with the
 165 c InsightII command file --- AL 1/1/95.
 166 c
 167             if (insight_cmd_out) then
 168             open (22,file=
 169      & 'insight_'//prefixp(:ilen(prefixp))//numm(:ilen(numm))//'.cmd',
 170      &   status='unknown')
 171 c
 172 c Write InsightII command file
 173 c
 174             cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_001'
 175             cff=ucase(cfname)
 176             write (22,'(5a)') 'get molecule pdb user ',
 177      &            cfname(:ilen(cfname)),'.pdb ',
 178      &            cff(:ilen(cff)),
 179      &            ' -heteroatom -reference_object'
 180             cfname=ucase(cfname)
 181             write (22,'(3a,i3,2(1h,,i3))') 'color molecule atoms ',
 182      &            cfname(:ilen(cfname)),' specified specification ',
 183      &            0,0,255
 184             write (22,'(2a)') 'display molecule only atoms heavy ',
 185      &            cfname(:ilen(cfname))
 186             cfname1=cfname
 187             deltae_max=energy(nconf(i,ncon_out))-emini
 188             do j=2,ncon_out
 189               deltae=energy(nconf(i,j))-emini
 190               icolor=aint(255*deltae/deltae_max)
 191               call numstr(j,mumm)
 192               cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_'
 193      &            //mumm(:ilen(mumm))
 194               cff=ucase(cfname)
 195               write (22,'(5a)') 'get molecule pdb user ',
 196      &            cfname(:ilen(cfname)),'.pdb ',
 197      &            cff(:ilen(cff)),' -heteroatom'
 198               cfname=ucase(cfname)
 199               do k=1,nct-nnt+1
 200                if (k.lt.10) then
 201                 write (22,'(3a,i1,3a,i1)')
 202      &        'superimpose -end_definition backbone "label mode" ',
 203      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 204      &         ':',k
 205                elseif (k.lt.100) then
 206                 write (22,'(3a,i2,3a,i2)')
 207      &        'superimpose -end_definition backbone "label mode" ',
 208      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 209      &         ':',k
 210                else
 211                 write (22,'(3a,i3,3a,i3)')
 212      &        'superimpose -end_definition backbone "label mode" ',
 213      &         cfname(:ilen(cfname)),':',k,' ',cfname1(:ilen(cfname1)),
 214      &         ':',k
 215                endif
 216
 217               enddo
 218               write (22,'(3a,i3,2(1h,,i3))') 'color molecule atoms ',
 219      &            cfname(:ilen(cfname)),' specified specification ',
 220      &            icolor,icolor,255-icolor
 221               write (22,'(2a)') 'display molecule only atoms heavy ',
 222      &            cfname(:ilen(cfname))
 223               close(22)
 224             enddo
 225             endif
 226
 227 c Write conformations of the family i to PDB files
 228             do j=1,ncon_out
 229               call numstr(j,mumm)
 230               icon=nconf(i,j)
 231               cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//'_'
 232      &             //mumm(:ilen(mumm))//exten
 233               do ii=1,nres
 234                 phi(ii)=phiall(ii,icon)
 235                 theta(ii)=thetall(ii,icon)
 236                 alph(ii)=alphall(ii,icon)
 237                 omeg(ii)=omall(ii,icon)
 238               enddo
 239               call chainbuild
 240               if (refstr.or.pdbref) rmsd=rmsnat(icon)
 241               open(ipdb,file=cfname,status='unknown',form='formatted')
 242               call pdbout(energy(icon),rmsd,titel)
 243               close(ipdb)
 244             enddo
 245           else
 246 c Produce only a single PDB file for the leading member of the family
 247             write (iout,*) 'Writing pdb file: icon=',icon
 248             if (from_cart .or. from_cx) then
 249
 250             do ii=1,2*nres
 251               do j=1,3
 252               c(j,ii)=allcart(j,ii,icon)
 253               enddo
 254             enddo
 255
 256             else
 257
 258             do ii=1,nres
 259               phi(ii)=phiall(ii,icon)
 260               theta(ii)=thetall(ii,icon)
 261               alph(ii)=alphall(ii,icon)
 262               omeg(ii)=omall(ii,icon)
 263             enddo
 264             call chainbuild
 265
 266             endif
 267
 268             cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//exten
 269             OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 270 c           print *,'Calling pdbout'
 271             if (refstr.or.pdbref) rmsd=rmsnat(icon)
 272             CALL PDBOUT(energy(icon),rmsd,titel)
 273             CLOSE(ipdb)
 274           endif
 275           I=I+1
 276           ICON=NCONF(I,1)
 277         ENDDO
 278       ENDIF
 279       IF (printmol2(icut).gt.0) THEN
 280 c Write out a number of conformations from each family in PDB format and
 281 c create InsightII command file for their displaying in different colors
 282         I=1
 283         ICON=NCONF(1,1)
 284         EMIN=ENERGY(ICON)
 285         DO WHILE(I.LE.NGR .AND. ENERGY(ICON)-EMIN.LE.ECUT)
 286 c          CALL NUMSTR(I,NUMM)
 287           write (NUMM,'(bz,i4.4)') i
 288           ncon_lim=min0(licz(i),printpdb(icut))
 289           write (iout,*) 'Writing mol2 file: icon=',icon
 290           do ii=1,nres
 291             phi(ii)=phiall(ii,icon)
 292             theta(ii)=thetall(ii,icon)
 293             alph(ii)=alphall(ii,icon)
 294             omeg(ii)=omall(ii,icon)
 295           enddo
 296           call chainbuild
 297           cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//extmol
 298           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 299 c         print *,'Calling pdbout'
 300           CALL MOL2OUT(energy(icon),'STRUCTURE'//numm)
 301           CLOSE(imol2)
 302           I=I+1
 303           ICON=NCONF(I,1)
 304         ENDDO
 305       ENDIF
 306   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 307   200 FORMAT (/'FAMILY ',I4,' CONTAINS ',I4,' CONFORMATION(S): ')
 308 c 300 FORMAT ( 8(I4,F6.1))
 309   300 FORMAT (5(I4,1pe12.4))
 310   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 311   500 FORMAT (8(2I4,2X))
 312   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 313       RETURN
 314       END