source/unres/src-5hdiag-tmp/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       include 'COMMON.LAGRANGE'
  13       character*50 tytul
  14       character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
  15       dimension ica(maxres)
  16       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
  17 cmodel      write (iunit,'(a5,i6)') 'MODEL',1
  18       if (nhfrag.gt.0) then
  19        do j=1,nhfrag
  20         iti=itype(hfrag(1,j))
  21         itj=itype(hfrag(2,j))
  22         if (j.lt.10) then
  23            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)')
  24      &           'HELIX',j,'H',j,
  25      &           restyp(iti),hfrag(1,j)-1,
  26      &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
  27         else
  28              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)')
  29      &           'HELIX',j,'H',j,
  30      &           restyp(iti),hfrag(1,j)-1,
  31      &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
  32         endif
  33        enddo
  34       endif
  35
  36       if (nbfrag.gt.0) then
  37
  38        do j=1,nbfrag
  39
  40         iti=itype(bfrag(1,j))
  41         itj=itype(bfrag(2,j)-1)
  42
  43         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)')
  44      &           'SHEET',1,'B',j,2,
  45      &           restyp(iti),bfrag(1,j)-1,
  46      &           restyp(itj),bfrag(2,j)-2,0
  47
  48         if (bfrag(3,j).gt.bfrag(4,j)) then
  49
  50          itk=itype(bfrag(3,j))
  51          itl=itype(bfrag(4,j)+1)
  52
  53          write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
  54      &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)')
  55      &           'SHEET',2,'B',j,2,
  56      &           restyp(itl),bfrag(4,j),
  57      &           restyp(itk),bfrag(3,j)-1,-1,
  58      &           "N",restyp(itk),bfrag(3,j)-1,
  59      &           "O",restyp(iti),bfrag(1,j)-1
  60
  61         else
  62
  63          itk=itype(bfrag(3,j))
  64          itl=itype(bfrag(4,j)-1)
  65
  66
  67         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
  68      &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)')
  69      &           'SHEET',2,'B',j,2,
  70      &           restyp(itk),bfrag(3,j)-1,
  71      &           restyp(itl),bfrag(4,j)-2,1,
  72      &           "N",restyp(itk),bfrag(3,j)-1,
  73      &           "O",restyp(iti),bfrag(1,j)-1
  74
  75
  76
  77         endif
  78
  79        enddo
  80       endif
  81
  82       if (nss.gt.0) then
  83         do i=1,nss
  84          if (dyn_ss) then
  85           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)')
  86      &         'SSBOND',i,'CYS',idssb(i)-nnt+1,
  87      &                    'CYS',jdssb(i)-nnt+1
  88          else
  89           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)')
  90      &         'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
  91      &                    'CYS',jhpb(i)-nnt+1-nres
  92          endif
  93         enddo
  94       endif
  95
  96       iatom=0
  97       ichain=1
  98       ires=0
  99       do i=nnt,nct
 100         iti=itype(i)
 101         if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then
 102           ichain=ichain+1
 103           ires=0
 104           write (iunit,'(a)') 'TER'
 105         else
 106         ires=ires+1
 107         iatom=iatom+1
 108         ica(i)=iatom
 109         if (iti.ne.ntyp1) then
 110         write (iunit,10) iatom,restyp(iti),chainid(ichain),
 111      &     ires,(c(j,i),j=1,3),vtot(i)
 112         if (iti.ne.10) then
 113           iatom=iatom+1
 114           write (iunit,20) iatom,restyp(iti),chainid(ichain),
 115      &      ires,(c(j,nres+i),j=1,3),
 116      &      vtot(i+nres)
 117          endif
 118         endif
 119         endif
 120       enddo
 121       write (iunit,'(a)') 'TER'
 122       do i=nnt,nct-1
 123         if (itype(i).eq.ntyp1) cycle
 124         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
 125           write (iunit,30) ica(i),ica(i+1)
 126         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
 127           write (iunit,30) ica(i),ica(i+1),ica(i)+1
 128         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
 129           write (iunit,30) ica(i),ica(i)+1
 130         endif
 131       enddo
 132       if (itype(nct).ne.10) then
 133         write (iunit,30) ica(nct),ica(nct)+1
 134       endif
 135       do i=1,nss
 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) || defined(CRAY)
 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       include 'COMMON.QRESTR'
 289       double precision time
 290 #if defined(AIX) || defined(PGI) || defined(CRAY)
 291       open(icart,file=cartname,position="append")
 292 #else
 293       open(icart,file=cartname,access="append")
 294 #endif
 295       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 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       include 'COMMON.QRESTR'
 332       double precision time
 333       integer iret,itmp
 334       real xcoord(3,maxres2+2),prec
 335
 336 #ifdef AIX
 337       call xdrfopen_(ixdrf,cartname, "a", iret)
 338       call xdrffloat_(ixdrf, real(time), iret)
 339       call xdrffloat_(ixdrf, real(potE), iret)
 340       call xdrffloat_(ixdrf, real(uconst), iret)
 341       call xdrffloat_(ixdrf, real(uconst_back), iret)
 342       call xdrffloat_(ixdrf, real(t_bath), iret)
 343       call xdrfint_(ixdrf, nss, iret)
 344       do j=1,nss
 345        if (dyn_ss) then
 346         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 347         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 348        else
 349         call xdrfint_(ixdrf, ihpb(j), iret)
 350         call xdrfint_(ixdrf, jhpb(j), iret)
 351        endif
 352       enddo
 353       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 354       do i=1,nfrag
 355         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 356       enddo
 357       do i=1,npair
 358         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 359       enddo
 360       do i=1,nfrag_back
 361         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 362         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 363         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 364       enddo
 365 #else
 366       call xdrfopen(ixdrf,cartname, "a", iret)
 367 c      write (iout,*) "Writing conformation: time",time," potE",potE,
 368 c     & " uconst",uconst," uconst_back",uconst_back," t_bath",t_bath,
 369 c     & " nss",nss
 370       call xdrffloat(ixdrf, real(time), iret)
 371       call xdrffloat(ixdrf, real(potE), iret)
 372       call xdrffloat(ixdrf, real(uconst), iret)
 373       call xdrffloat(ixdrf, real(uconst_back), iret)
 374       call xdrffloat(ixdrf, real(t_bath), iret)
 375       call xdrfint(ixdrf, nss, iret)
 376       do j=1,nss
 377        if (dyn_ss) then
 378         call xdrfint(ixdrf, idssb(j)+nres, iret)
 379         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 380        else
 381         call xdrfint(ixdrf, ihpb(j), iret)
 382         call xdrfint(ixdrf, jhpb(j), iret)
 383        endif
 384       enddo
 385       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 386       do i=1,nfrag
 387         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 388       enddo
 389       do i=1,npair
 390         call xdrffloat(ixdrf, real(qpair(i)), iret)
 391       enddo
 392       do i=1,nfrag_back
 393         call xdrffloat(ixdrf, real(utheta(i)), iret)
 394         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 395         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 396       enddo
 397 #endif
 398       prec=10000.0
 399       do i=1,nres
 400        do j=1,3
 401         xcoord(j,i)=c(j,i)
 402        enddo
 403       enddo
 404       do i=nnt,nct
 405        do j=1,3
 406         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 407        enddo
 408       enddo
 409
 410       itmp=nres+nct-nnt+1
 411 #ifdef AIX
 412       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 413       call xdrfclose_(ixdrf, iret)
 414 #else
 415       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 416       call xdrfclose(ixdrf, iret)
 417 #endif
 418       return
 419       end
 420 #endif
 421 c-----------------------------------------------------------------
 422       subroutine statout(itime)
 423       implicit real*8 (a-h,o-z)
 424       include 'DIMENSIONS'
 425       include 'COMMON.CONTROL'
 426       include 'COMMON.CHAIN'
 427       include 'COMMON.INTERACT'
 428       include 'COMMON.NAMES'
 429       include 'COMMON.IOUNITS'
 430       include 'COMMON.HEADER'
 431       include 'COMMON.SBRIDGE'
 432       include 'COMMON.DISTFIT'
 433       include 'COMMON.MD'
 434       include 'COMMON.QRESTR'
 435       include 'COMMON.REMD'
 436       include 'COMMON.SETUP'
 437       integer itime
 438       double precision energia(0:n_ene)
 439       double precision gyrate
 440       external gyrate
 441       common /gucio/ cm
 442       character*256 line1,line2
 443       character*4 format1,format2
 444       character*30 format
 445 #ifdef AIX
 446       if(itime.eq.0) then
 447        open(istat,file=statname,position="append")
 448       endif
 449 #else
 450 #if defined(PGI) || defined(CRAY)
 451       open(istat,file=statname,position="append")
 452 #else
 453       open(istat,file=statname,access="append")
 454 #endif
 455 #endif
 456        if (AFMlog.gt.0) then
 457        if (refstr) then
 458          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 459           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')
 460      &          itime,totT,EK,potE,totE,
 461      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
 462      &          potEcomp(23),me
 463           format1="a133"
 464          else
 465 C          print *,'A CHUJ',potEcomp(23)
 466           write (line1,'(i10,f15.2,7f12.3,i5,$)')
 467      &           itime,totT,EK,potE,totE,
 468      &           kinetic_T,t_bath,gyrate(),
 469      &           potEcomp(23),me
 470           format1="a114"
 471         endif
 472        else if (selfguide.gt.0) then
 473        distance=0.0
 474        do j=1,3
 475        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
 476        enddo
 477        distance=dsqrt(distance)
 478        if (refstr) then
 479          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 480           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2,
 481      &    f9.3,i5,$)')
 482      &          itime,totT,EK,potE,totE,
 483      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
 484      &          distance,potEcomp(23),me
 485           format1="a133"
 486 C          print *,"CHUJOWO"
 487          else
 488 C          print *,'A CHUJ',potEcomp(23)
 489           write (line1,'(i10,f15.2,8f12.3,i5,$)')
 490      &           itime,totT,EK,potE,totE,
 491      &           kinetic_T,t_bath,gyrate(),
 492      &           distance,potEcomp(23),me
 493           format1="a114"
 494         endif
 495        else
 496        if (refstr) then
 497          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 498           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')
 499      &          itime,totT,EK,potE,totE,
 500      &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 501           format1="a133"
 502         else
 503           write (line1,'(i10,f15.2,7f12.3,i5,$)')
 504      &           itime,totT,EK,potE,totE,
 505      &           amax,kinetic_T,t_bath,gyrate(),me
 506           format1="a114"
 507         endif
 508         endif
 509         if(usampl.and.totT.gt.eq_time) then
 510            if (loc_qlike) then
 511            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
 512      &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
 513      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back),
 514      &      ((qloc(j,i),j=1,3),i=1,nfrag_back)
 515            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
 516      &             +42*nfrag_back
 517            else
 518            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
 519      &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
 520      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 521            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
 522      &             +21*nfrag_back
 523            endif
 524         else
 525            format2="a001"
 526            line2=' '
 527         endif
 528         if (print_compon) then
 529           if(itime.eq.0) then
 530            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
 531      &                                                     ",31a12)"
 532            write (istat,format) "#","",
 533      &      (ename(print_order(i)),i=1,nprint_ene)
 534           endif
 535           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
 536      &                                                     ",31f12.3)"
 537           write (istat,format) line1,line2,
 538      &      (potEcomp(print_order(i)),i=1,nprint_ene)
 539         else
 540           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 541           write (istat,format) line1,line2
 542         endif
 543 #if defined(AIX)
 544         call flush(istat)
 545 #else
 546         close(istat)
 547 #endif
 548        return
 549       end
 550 c---------------------------------------------------------------
 551       double precision function gyrate()
 552       implicit real*8 (a-h,o-z)
 553       include 'DIMENSIONS'
 554       include 'COMMON.INTERACT'
 555       include 'COMMON.CHAIN'
 556       double precision cen(3),rg
 557
 558       do j=1,3
 559        cen(j)=0.0d0
 560       enddo
 561
 562       ii=0
 563       do i=nnt,nct
 564         if (itype(i).eq.ntyp1) cycle
 565         ii=ii+1
 566         do j=1,3
 567           cen(j)=cen(j)+c(j,i)
 568         enddo
 569       enddo
 570       do j=1,3
 571         cen(j)=cen(j)/dble(ii)
 572       enddo
 573       rg = 0.0d0
 574       do i = nnt, nct
 575         if (itype(i).eq.ntyp1) cycle
 576         do j=1,3
 577          rg = rg + (c(j,i)-cen(j))**2
 578         enddo
 579       end do
 580       gyrate = dsqrt(rg/dble(ii))
 581       return
 582       end