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