94fd7d747970e7f17d4f59908e31c2c5553a200b
[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*(*) 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       write (iout,*) "cartout: cartname ",cartname
272 #if defined(AIX) || defined(PGI)
273       open(icart,file=cartname,position="append")
274 #else
275       open(icart,file=cartname,access="append")
276 #endif
277       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
278       if (dyn_ss) then
279        write (icart,'(i4,$)')
280      &   nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)       
281       else
282        write (icart,'(i4,$)')
283      &   nss,(ihpb(j),jhpb(j),j=1,nss)
284        endif
285        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,
286      & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),
287      & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
288       write (icart,'(8f10.5)')
289      & ((c(k,j),k=1,3),j=1,nres),
290      & ((c(k,j+nres),k=1,3),j=nnt,nct)
291       close(icart)
292       return
293       end
294 c-----------------------------------------------------------------
295 #ifndef NOXDR
296       subroutine cartout(time)
297       implicit real*8 (a-h,o-z)
298       include 'DIMENSIONS'
299 #ifdef MPI
300       include 'mpif.h'
301       include 'COMMON.SETUP'
302 #else
303       parameter (me=0)
304 #endif
305       include 'COMMON.CHAIN'
306       include 'COMMON.INTERACT'
307       include 'COMMON.NAMES'
308       include 'COMMON.IOUNITS'
309       include 'COMMON.HEADER'
310       include 'COMMON.SBRIDGE'
311       include 'COMMON.DISTFIT'
312       include 'COMMON.MD'
313       double precision time
314       integer iret,itmp
315       real xcoord(3,maxres2+2),prec
316
317 #ifdef AIX
318       call xdrfopen_(ixdrf,cartname, "a", iret)
319       call xdrffloat_(ixdrf, real(time), iret)
320       call xdrffloat_(ixdrf, real(potE), iret)
321       call xdrffloat_(ixdrf, real(uconst), iret)
322       call xdrffloat_(ixdrf, real(uconst_back), iret)
323       call xdrffloat_(ixdrf, real(t_bath), iret)
324       call xdrfint_(ixdrf, nss, iret) 
325       do j=1,nss
326        if (dyn_ss) then
327         call xdrfint_(ixdrf, idssb(j)+nres, iret)
328         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
329        else
330         call xdrfint_(ixdrf, ihpb(j), iret)
331         call xdrfint_(ixdrf, jhpb(j), iret)
332        endif
333       enddo
334       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
335       do i=1,nfrag
336         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
337       enddo
338       do i=1,npair
339         call xdrffloat_(ixdrf, real(qpair(i)), iret)
340       enddo
341       do i=1,nfrag_back
342         call xdrffloat_(ixdrf, real(utheta(i)), iret)
343         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
344         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
345       enddo
346 #else
347       call xdrfopen(ixdrf,cartname, "a", iret)
348       call xdrffloat(ixdrf, real(time), iret)
349       call xdrffloat(ixdrf, real(potE), iret)
350       call xdrffloat(ixdrf, real(uconst), iret)
351       call xdrffloat(ixdrf, real(uconst_back), iret)
352       call xdrffloat(ixdrf, real(t_bath), iret)
353       call xdrfint(ixdrf, nss, iret) 
354       do j=1,nss
355        if (dyn_ss) then
356         call xdrfint(ixdrf, idssb(j)+nres, iret)
357         call xdrfint(ixdrf, jdssb(j)+nres, iret)
358        else
359         call xdrfint(ixdrf, ihpb(j), iret)
360         call xdrfint(ixdrf, jhpb(j), iret)
361        endif
362       enddo
363       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
364       do i=1,nfrag
365         call xdrffloat(ixdrf, real(qfrag(i)), iret)
366       enddo
367       do i=1,npair
368         call xdrffloat(ixdrf, real(qpair(i)), iret)
369       enddo
370       do i=1,nfrag_back
371         call xdrffloat(ixdrf, real(utheta(i)), iret)
372         call xdrffloat(ixdrf, real(ugamma(i)), iret)
373         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
374       enddo
375 #endif
376       prec=10000.0
377       do i=1,nres
378        do j=1,3
379         xcoord(j,i)=c(j,i)
380        enddo
381       enddo
382       do i=nnt,nct
383        do j=1,3
384         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
385        enddo
386       enddo
387
388       itmp=nres+nct-nnt+1
389 #ifdef AIX
390       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
391       call xdrfclose_(ixdrf, iret)
392 #else
393       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
394       call xdrfclose(ixdrf, iret)
395 #endif
396       return
397       end
398 #endif
399 c-----------------------------------------------------------------
400       subroutine statout(itime)
401       implicit real*8 (a-h,o-z)
402       include 'DIMENSIONS'
403       include 'COMMON.CONTROL'
404       include 'COMMON.CHAIN'
405       include 'COMMON.INTERACT'
406       include 'COMMON.NAMES'
407       include 'COMMON.IOUNITS'
408       include 'COMMON.HEADER'
409       include 'COMMON.SBRIDGE'
410       include 'COMMON.DISTFIT'
411       include 'COMMON.MD'
412       include 'COMMON.REMD'
413       include 'COMMON.SETUP'
414       integer itime
415       double precision energia(0:n_ene)
416       double precision gyrate
417       external gyrate
418       common /gucio/ cm
419       character*256 line1,line2
420       character*4 format1,format2
421       character*30 format
422 #ifdef AIX
423       if(itime.eq.0) then
424        open(istat,file=statname,position="append")
425       endif
426 #else
427 #ifdef PGI
428       open(istat,file=statname,position="append")
429 #else
430       open(istat,file=statname,access="append")
431 #endif
432 #endif
433        if (refstr) then
434          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
435         if(tnp .or. tnp1 .or. tnh) then
436         write (line1,'(i10,f15.2,3f12.3,f12.6,f7.2,4f6.3,3f12.3,i5,$)')
437      &          itime,totT,EK,potE,totE,hhh,
438      &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
439           format1="a145"
440         else
441           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')
442      &          itime,totT,EK,potE,totE,
443      &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
444           format1="a133"
445         endif
446        else
447         if(tnp .or. tnp1 .or. tnh) then
448           write (line1,'(i10,f15.2,7f12.3,f12.6,i5,$)')
449      &           itime,totT,EK,potE,totE,hhh,
450      &           amax,kinetic_T,t_bath,gyrate(),me
451           format1="a126"
452         else
453           write (line1,'(i10,f15.2,7f12.3,i5,$)')
454      &           itime,totT,EK,potE,totE,
455      &           amax,kinetic_T,t_bath,gyrate(),me
456           format1="a114"
457         endif
458        endif
459         if(usampl.and.totT.gt.eq_time) then
460            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
461      &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
462      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
463            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
464      &             +21*nfrag_back
465         elseif(hremd.gt.0.or.homol_nset.gt.1) then
466            write(line2,'(i5)') iset
467            format2="a005"
468         else
469            format2="a001"
470            line2=' '
471         endif
472         if (print_compon) then
473           if(itime.eq.0) then
474            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
475      &                                                     ",100a12)"
476            write (istat,format) "#","",
477      &      (ename(print_order(i)),i=1,nprint_ene)
478           endif
479           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
480      &                                                     ",100f12.3)"
481           write (istat,format) line1,line2,
482      &      (potEcomp(print_order(i)),i=1,nprint_ene)
483         else
484           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
485           write (istat,format) line1,line2
486         endif
487 #if defined(AIX)
488         call flush(istat)
489 #else
490         close(istat)
491 #endif
492        return
493       end
494 c---------------------------------------------------------------                      
495       double precision function gyrate()
496       implicit real*8 (a-h,o-z)
497       include 'DIMENSIONS'
498       include 'COMMON.INTERACT'
499       include 'COMMON.CHAIN'
500       double precision cen(3),rg
501
502       do j=1,3
503        cen(j)=0.0d0
504       enddo
505
506       do i=nnt,nct
507           do j=1,3
508             cen(j)=cen(j)+c(j,i)
509           enddo
510       enddo
511       do j=1,3
512             cen(j)=cen(j)/dble(nct-nnt+1)
513       enddo
514       rg = 0.0d0
515       do i = nnt, nct
516         do j=1,3
517          rg = rg + (c(j,i)-cen(j))**2 
518         enddo
519       end do
520       gyrate = sqrt(rg/dble(nct-nnt+1))
521       return
522       end
523