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