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