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