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