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