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