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