wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / geomout.F
1       subroutine pdbout(etot,tytul,iunit)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.CHAIN'
5       include 'COMMON.INTERACT'
6       include 'COMMON.NAMES'
7       include 'COMMON.IOUNITS'
8       include 'COMMON.HEADER'
9       include 'COMMON.SBRIDGE'
10       include 'COMMON.DISTFIT'
11       include 'COMMON.MD'
12       character*50 tytul
13       character*1 chainid(26) /'A','B','C','D','E','F','G','H','I','J',
14      &  'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
15
16       dimension ica(maxres)
17       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
18 cmodel      write (iunit,'(a5,i6)') 'MODEL',1
19       if (nhfrag.gt.0) then
20        do j=1,nhfrag
21         iti=itype(hfrag(1,j))
22         itj=itype(hfrag(2,j))
23         if (j.lt.10) then
24            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') 
25      &           'HELIX',j,'H',j,
26      &           restyp(iti),hfrag(1,j)-1,
27      &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
28         else
29              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') 
30      &           'HELIX',j,'H',j,
31      &           restyp(iti),hfrag(1,j)-1,
32      &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
33         endif
34        enddo
35       endif
36
37       if (nbfrag.gt.0) then
38
39        do j=1,nbfrag
40
41         iti=itype(bfrag(1,j))
42         itj=itype(bfrag(2,j)-1)
43
44         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') 
45      &           'SHEET',1,'B',j,2,
46      &           restyp(iti),bfrag(1,j)-1,
47      &           restyp(itj),bfrag(2,j)-2,0
48
49         if (bfrag(3,j).gt.bfrag(4,j)) then
50
51          itk=itype(bfrag(3,j))
52          itl=itype(bfrag(4,j)+1)
53
54          write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
55      &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') 
56      &           'SHEET',2,'B',j,2,
57      &           restyp(itl),bfrag(4,j),
58      &           restyp(itk),bfrag(3,j)-1,-1,
59      &           "N",restyp(itk),bfrag(3,j)-1,
60      &           "O",restyp(iti),bfrag(1,j)-1
61
62         else
63
64          itk=itype(bfrag(3,j))
65          itl=itype(bfrag(4,j)-1)
66
67
68         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
69      &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') 
70      &           'SHEET',2,'B',j,2,
71      &           restyp(itk),bfrag(3,j)-1,
72      &           restyp(itl),bfrag(4,j)-2,1,
73      &           "N",restyp(itk),bfrag(3,j)-1,
74      &           "O",restyp(iti),bfrag(1,j)-1
75
76
77
78         endif
79          
80        enddo
81       endif 
82
83       if (nss.gt.0) then
84         do i=1,nss
85          if (dyn_ss) then
86           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
87      &         'SSBOND',i,'CYS',idssb(i)-nnt+1,
88      &                    'CYS',jdssb(i)-nnt+1
89          else
90           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
91      &         'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
92      &                    'CYS',jhpb(i)-nnt+1-nres
93          endif
94         enddo
95       endif
96       
97       iatom=0
98       ichain=1
99       ires=0
100       do i=nnt,nct
101         iti=itype(i)
102         if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then
103           ichain=ichain+1
104           ires=0
105           write (iunit,'(a)') 'TER'
106         else
107         ires=ires+1
108         iatom=iatom+1
109         ica(i)=iatom
110         if (iti.ne.ntyp1) then
111         write (iunit,10) iatom,restyp(iti),chainid(1+mod(ichain/2,26)),
112      &     ires,(c(j,i),j=1,3),vtot(i)
113         if (iti.ne.10) then
114           iatom=iatom+1
115         write (iunit,20) iatom,restyp(iti),chainid(1+mod(ichain/2,26)),
116      &      ires,(c(j,nres+i),j=1,3),
117      &      vtot(i+nres)
118         endif
119         endif
120         endif
121       enddo
122       write (iunit,'(a)') 'TER'
123       do i=nnt,nct-1
124         if (itype(i).eq.ntyp1) cycle
125         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
126           write (iunit,30) ica(i),ica(i+1)
127         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
128           write (iunit,30) ica(i),ica(i+1),ica(i)+1
129         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
130           write (iunit,30) ica(i),ica(i)+1
131         endif
132       enddo
133       if (itype(nct).ne.10) then
134         write (iunit,30) ica(nct),ica(nct)+1
135       endif
136       do i=1,nss
137        if (dyn_ss) then
138         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
139        else
140         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
141        endif
142       enddo
143       write (iunit,'(a6)') 'ENDMDL'     
144   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
145   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
146   30  FORMAT ('CONECT',8I5)
147       return
148       end
149 c------------------------------------------------------------------------------
150       subroutine MOL2out(etot,tytul)
151 C Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2 
152 C format.
153       implicit real*8 (a-h,o-z)
154       include 'DIMENSIONS'
155       include 'COMMON.CHAIN'
156       include 'COMMON.INTERACT'
157       include 'COMMON.NAMES'
158       include 'COMMON.IOUNITS'
159       include 'COMMON.HEADER'
160       include 'COMMON.SBRIDGE'
161       character*50 tytul,fd
162       character*3 zahl
163       character*6 res_num,pom,ucase
164 #ifdef AIX
165       call fdate_(fd)
166 #elif (defined CRAY)
167       call date(fd)
168 #else
169       call fdate(fd)
170 #endif
171       write (imol2,'(a)') '#'
172       write (imol2,'(a)') 
173      & '#         Creating user name:           unres'
174       write (imol2,'(2a)') '#         Creation time:                ',
175      & fd
176       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
177       write (imol2,'(a)') tytul
178       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
179       write (imol2,'(a)') 'SMALL'
180       write (imol2,'(a)') 'USER_CHARGES'
181       write (imol2,'(a)') '\@<TRIPOS>ATOM' 
182       do i=nnt,nct
183         write (zahl,'(i3)') i
184         pom=ucase(restyp(itype(i)))
185         res_num = pom(:3)//zahl(2:)
186         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
187       enddo
188       write (imol2,'(a)') '\@<TRIPOS>BOND'
189       do i=nnt,nct-1
190         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
191       enddo
192       do i=1,nss
193         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
194       enddo
195       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
196       do i=nnt,nct
197         write (zahl,'(i3)') i
198         pom = ucase(restyp(itype(i)))
199         res_num = pom(:3)//zahl(2:)
200         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
201       enddo
202   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
203   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
204       return
205       end
206 c------------------------------------------------------------------------
207       subroutine intout
208       implicit real*8 (a-h,o-z)
209       include 'DIMENSIONS'
210       include 'COMMON.IOUNITS'
211       include 'COMMON.CHAIN'
212       include 'COMMON.VAR'
213       include 'COMMON.LOCAL'
214       include 'COMMON.INTERACT'
215       include 'COMMON.NAMES'
216       include 'COMMON.GEO'
217       include 'COMMON.TORSION'
218       write (iout,'(/a)') 'Geometry of the virtual chain.'
219       write (iout,'(7a)') '  Res  ','         d','     Theta',
220      & '     Gamma','       Dsc','     Alpha','      Beta '
221       do i=1,nres
222         iti=itype(i)
223         write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),
224      &     rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),
225      &     rad2deg*omeg(i)
226       enddo
227       return
228       end
229 c---------------------------------------------------------------------------
230       subroutine briefout(it,ener)
231       implicit real*8 (a-h,o-z)
232       include 'DIMENSIONS'
233       include 'COMMON.IOUNITS'
234       include 'COMMON.CHAIN'
235       include 'COMMON.VAR'
236       include 'COMMON.LOCAL'
237       include 'COMMON.INTERACT'
238       include 'COMMON.NAMES'
239       include 'COMMON.GEO'
240       include 'COMMON.SBRIDGE'
241 c     print '(a,i5)',intname,igeom
242 #if defined(AIX) || defined(PGI)
243       open (igeom,file=intname,position='append')
244 #else
245       open (igeom,file=intname,access='append')
246 #endif
247       IF (NSS.LE.9) THEN
248         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
249       ELSE
250         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
251         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
252       ENDIF
253 c     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
254       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
255       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
256 c     if (nvar.gt.nphi+ntheta) then
257         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
258         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
259 c     endif
260       close(igeom)
261   180 format (I5,F12.3,I2,9(1X,2I3))
262   190 format (3X,11(1X,2I3))
263   200 format (8F10.4)
264       return
265       end
266 #ifdef WINIFL
267       subroutine fdate(fd)
268       character*32 fd
269       write(fd,'(32x)')
270       return
271       end
272 #endif
273 c----------------------------------------------------------------
274 #ifdef NOXDR
275       subroutine cartout(time)
276 #else
277       subroutine cartoutx(time)
278 #endif
279       implicit real*8 (a-h,o-z)
280       include 'DIMENSIONS'
281       include 'COMMON.CHAIN'
282       include 'COMMON.INTERACT'
283       include 'COMMON.NAMES'
284       include 'COMMON.IOUNITS'
285       include 'COMMON.HEADER'
286       include 'COMMON.SBRIDGE'
287       include 'COMMON.DISTFIT'
288       include 'COMMON.MD'
289       double precision time
290       write (iout,*) "cartout: cartname ",cartname
291 #if defined(AIX) || defined(PGI)
292       open(icart,file=cartname,position="append")
293 #else
294       open(icart,file=cartname,access="append")
295 #endif
296       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
297       if (dyn_ss) then
298        write (icart,'(i4,$)')
299      &   nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)       
300       else
301        write (icart,'(i4,$)')
302      &   nss,(ihpb(j),jhpb(j),j=1,nss)
303        endif
304        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,
305      & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),
306      & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
307       write (icart,'(8f10.5)')
308      & ((c(k,j),k=1,3),j=1,nres),
309      & ((c(k,j+nres),k=1,3),j=nnt,nct)
310       close(icart)
311       return
312       end
313 c-----------------------------------------------------------------
314 #ifndef NOXDR
315       subroutine cartout(time)
316       implicit real*8 (a-h,o-z)
317       include 'DIMENSIONS'
318 #ifdef MPI
319       include 'mpif.h'
320       include 'COMMON.SETUP'
321 #else
322       parameter (me=0)
323 #endif
324       include 'COMMON.CHAIN'
325       include 'COMMON.INTERACT'
326       include 'COMMON.NAMES'
327       include 'COMMON.IOUNITS'
328       include 'COMMON.HEADER'
329       include 'COMMON.SBRIDGE'
330       include 'COMMON.DISTFIT'
331       include 'COMMON.MD'
332       double precision time
333       integer iret,itmp
334       real xcoord(3,maxres2+2),prec
335
336 #ifdef AIX
337       call xdrfopen_(ixdrf,cartname, "a", iret)
338       call xdrffloat_(ixdrf, real(time), iret)
339       call xdrffloat_(ixdrf, real(potE), iret)
340       call xdrffloat_(ixdrf, real(uconst), iret)
341       call xdrffloat_(ixdrf, real(uconst_back), iret)
342       call xdrffloat_(ixdrf, real(t_bath), iret)
343       call xdrfint_(ixdrf, nss, iret) 
344       do j=1,nss
345        if (dyn_ss) then
346         call xdrfint_(ixdrf, idssb(j)+nres, iret)
347         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
348        else
349         call xdrfint_(ixdrf, ihpb(j), iret)
350         call xdrfint_(ixdrf, jhpb(j), iret)
351        endif
352       enddo
353       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
354       do i=1,nfrag
355         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
356       enddo
357       do i=1,npair
358         call xdrffloat_(ixdrf, real(qpair(i)), iret)
359       enddo
360       do i=1,nfrag_back
361         call xdrffloat_(ixdrf, real(utheta(i)), iret)
362         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
363         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
364       enddo
365 #else
366       call xdrfopen(ixdrf,cartname, "a", iret)
367       call xdrffloat(ixdrf, real(time), iret)
368       call xdrffloat(ixdrf, real(potE), iret)
369       call xdrffloat(ixdrf, real(uconst), iret)
370       call xdrffloat(ixdrf, real(uconst_back), iret)
371       call xdrffloat(ixdrf, real(t_bath), iret)
372       call xdrfint(ixdrf, nss, iret) 
373       do j=1,nss
374        if (dyn_ss) then
375         call xdrfint(ixdrf, idssb(j)+nres, iret)
376         call xdrfint(ixdrf, jdssb(j)+nres, iret)
377        else
378         call xdrfint(ixdrf, ihpb(j), iret)
379         call xdrfint(ixdrf, jhpb(j), iret)
380        endif
381       enddo
382       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
383       do i=1,nfrag
384         call xdrffloat(ixdrf, real(qfrag(i)), iret)
385       enddo
386       do i=1,npair
387         call xdrffloat(ixdrf, real(qpair(i)), iret)
388       enddo
389       do i=1,nfrag_back
390         call xdrffloat(ixdrf, real(utheta(i)), iret)
391         call xdrffloat(ixdrf, real(ugamma(i)), iret)
392         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
393       enddo
394 #endif
395       prec=10000.0
396       do i=1,nres
397        do j=1,3
398         xcoord(j,i)=c(j,i)
399        enddo
400       enddo
401       do i=nnt,nct
402        do j=1,3
403         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
404        enddo
405       enddo
406
407       itmp=nres+nct-nnt+1
408 #ifdef AIX
409       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
410       call xdrfclose_(ixdrf, iret)
411 #else
412       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
413       call xdrfclose(ixdrf, iret)
414 #endif
415       return
416       end
417 #endif
418 c-----------------------------------------------------------------
419       subroutine statout(itime)
420       implicit real*8 (a-h,o-z)
421       include 'DIMENSIONS'
422       include 'COMMON.CONTROL'
423       include 'COMMON.CHAIN'
424       include 'COMMON.INTERACT'
425       include 'COMMON.NAMES'
426       include 'COMMON.IOUNITS'
427       include 'COMMON.HEADER'
428       include 'COMMON.SBRIDGE'
429       include 'COMMON.DISTFIT'
430       include 'COMMON.MD'
431       include 'COMMON.REMD'
432       include 'COMMON.SETUP'
433       integer itime
434       double precision energia(0:n_ene)
435       double precision gyrate
436       external gyrate
437       common /gucio/ cm
438       character*256 line1,line2
439       character*4 format1,format2
440       character*30 format
441 #ifdef AIX
442       if(itime.eq.0) then
443        open(istat,file=statname,position="append")
444       endif
445 #else
446 #ifdef PGI
447       open(istat,file=statname,position="append")
448 #else
449       open(istat,file=statname,access="append")
450 #endif
451 #endif
452        if (AFMlog.gt.0) then
453        if (refstr) then
454          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
455           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')
456      &          itime,totT,EK,potE,totE,
457      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
458      &          potEcomp(23),me
459           format1="a133"
460        else
461 C          print *,'A CHUJ',potEcomp(23)
462           write (line1,'(i10,f15.2,7f12.3,i5,$)')
463      &           itime,totT,EK,potE,totE,
464      &           kinetic_T,t_bath,gyrate(),
465      &           potEcomp(23),me
466           format1="a114"
467         endif
468        else if (selfguide.gt.0) then
469        distance=0.0
470        do j=1,3
471        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
472        enddo
473        distance=dsqrt(distance)
474        if (refstr) then
475          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
476           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2,
477      &    f9.2,i5,$)')
478      &          itime,totT,EK,potE,totE,
479      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
480      &          distance,potEcomp(23),me
481           format1="a133"
482 C          print *,"CHUJOWO"
483          else
484 C          print *,'A CHUJ',potEcomp(23)
485           write (line1,'(i10,f15.2,8f12.3,i5,$)')
486      &           itime,totT,EK,potE,totE,
487      &           kinetic_T,t_bath,gyrate(),
488      &           distance,potEcomp(23),me
489           format1="a114"
490         endif
491        else
492         if (refstr) then
493          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
494           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')
495      &          itime,totT,EK,potE,totE,
496      &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
497           format1="a133"
498         else
499           write (line1,'(i10,f15.2,7f12.3,i5,$)')
500      &           itime,totT,EK,potE,totE,
501      &           amax,kinetic_T,t_bath,gyrate(),me
502           format1="a114"
503         endif
504        endif
505         if(usampl.and.totT.gt.eq_time) then
506            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
507      &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
508      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
509            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
510      &             +21*nfrag_back
511         elseif(hremd.gt.0.or.homol_nset.gt.1) then
512            write(line2,'(i5)') iset
513            format2="a005"
514         else
515            format2="a001"
516            line2=' '
517         endif
518         if (print_compon) then
519 #ifdef DEBUG
520           write (iout,*) "itime",itime," temperature",t_bath,
521      &     " potential energy",potE,potEcomp(0)
522           call enerprint(potEcomp)
523 #endif
524           if(itime.eq.0) then
525            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
526      &                                                     ",100a12)"
527            write (istat,format) "#","",
528      &      (ename(print_order(i)),i=1,nprint_ene)
529           endif
530           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
531      &                                                     ",100f12.3)"
532           write (istat,format) line1,line2,
533      &      (potEcomp(print_order(i)),i=1,nprint_ene)
534         else
535           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
536           write (istat,format) line1,line2
537         endif
538 #if defined(AIX)
539         call flush(istat)
540 #else
541         close(istat)
542 #endif
543        return
544       end
545 c---------------------------------------------------------------                      
546       double precision function gyrate()
547       implicit real*8 (a-h,o-z)
548       include 'DIMENSIONS'
549       include 'COMMON.INTERACT'
550       include 'COMMON.CHAIN'
551       double precision cen(3),rg
552
553       do j=1,3
554        cen(j)=0.0d0
555       enddo
556
557       do i=nnt,nct
558           do j=1,3
559             cen(j)=cen(j)+c(j,i)
560           enddo
561       enddo
562       do j=1,3
563             cen(j)=cen(j)/dble(nct-nnt+1)
564       enddo
565       rg = 0.0d0
566       do i = nnt, nct
567         do j=1,3
568          rg = rg + (c(j,i)-cen(j))**2 
569         enddo
570       end do
571       gyrate = sqrt(rg/dble(nct-nnt+1))
572       return
573       end
574