force legacy standard for gfortran
[unres.git] / 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