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*32 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      & '       Phi','       Dsc','     Alpha','      Omega'
 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 #if defined(AIX) || defined(PGI)
 289       open(icart,file=cartname,position="append")
 290 #else
 291       open(icart,file=cartname,access="append")
 292 #endif
 293       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 294       if (dyn_ss) then
 295        write (icart,'(i4,$)')
 296      &   nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
 297       else
 298        write (icart,'(i4,$)')
 299      &   nss,(ihpb(j),jhpb(j),j=1,nss)
 300        endif
 301        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,
 302      & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),
 303      & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 304       write (icart,'(8f10.5)')
 305      & ((c(k,j),k=1,3),j=1,nres),
 306      & ((c(k,j+nres),k=1,3),j=nnt,nct)
 307       close(icart)
 308       return
 309       end
 310 c-----------------------------------------------------------------
 311 #ifndef NOXDR
 312       subroutine cartout(time)
 313       implicit real*8 (a-h,o-z)
 314       include 'DIMENSIONS'
 315 #ifdef MPI
 316       include 'mpif.h'
 317       include 'COMMON.SETUP'
 318 #else
 319       parameter (me=0)
 320 #endif
 321       include 'COMMON.CHAIN'
 322       include 'COMMON.INTERACT'
 323       include 'COMMON.NAMES'
 324       include 'COMMON.IOUNITS'
 325       include 'COMMON.HEADER'
 326       include 'COMMON.SBRIDGE'
 327       include 'COMMON.DISTFIT'
 328       include 'COMMON.MD'
 329       double precision time
 330       integer iret,itmp
 331       real xcoord(3,maxres2+2),prec
 332
 333 #ifdef AIX
 334       call xdrfopen_(ixdrf,cartname, "a", iret)
 335       call xdrffloat_(ixdrf, real(time), iret)
 336       call xdrffloat_(ixdrf, real(potE), iret)
 337       call xdrffloat_(ixdrf, real(uconst), iret)
 338       call xdrffloat_(ixdrf, real(uconst_back), iret)
 339       call xdrffloat_(ixdrf, real(t_bath), iret)
 340       call xdrfint_(ixdrf, nss, iret)
 341       do j=1,nss
 342        if (dyn_ss) then
 343         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 344         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 345        else
 346         call xdrfint_(ixdrf, ihpb(j), iret)
 347         call xdrfint_(ixdrf, jhpb(j), iret)
 348        endif
 349       enddo
 350       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 351       do i=1,nfrag
 352         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 353       enddo
 354       do i=1,npair
 355         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 356       enddo
 357       do i=1,nfrag_back
 358         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 359         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 360         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 361       enddo
 362 #else
 363       call xdrfopen(ixdrf,cartname, "a", iret)
 364       call xdrffloat(ixdrf, real(time), iret)
 365       call xdrffloat(ixdrf, real(potE), iret)
 366       call xdrffloat(ixdrf, real(uconst), iret)
 367       call xdrffloat(ixdrf, real(uconst_back), iret)
 368       call xdrffloat(ixdrf, real(t_bath), iret)
 369       call xdrfint(ixdrf, nss, iret)
 370       do j=1,nss
 371        if (dyn_ss) then
 372         call xdrfint(ixdrf, idssb(j)+nres, iret)
 373         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 374        else
 375         call xdrfint(ixdrf, ihpb(j), iret)
 376         call xdrfint(ixdrf, jhpb(j), iret)
 377        endif
 378       enddo
 379       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 380       do i=1,nfrag
 381         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 382       enddo
 383       do i=1,npair
 384         call xdrffloat(ixdrf, real(qpair(i)), iret)
 385       enddo
 386       do i=1,nfrag_back
 387         call xdrffloat(ixdrf, real(utheta(i)), iret)
 388         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 389         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 390       enddo
 391 #endif
 392       prec=10000.0
 393       do i=1,nres
 394        do j=1,3
 395         xcoord(j,i)=c(j,i)
 396        enddo
 397       enddo
 398       do i=nnt,nct
 399        do j=1,3
 400         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 401        enddo
 402       enddo
 403
 404       itmp=nres+nct-nnt+1
 405 #ifdef AIX
 406       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 407       call xdrfclose_(ixdrf, iret)
 408 #else
 409       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 410       call xdrfclose(ixdrf, iret)
 411 #endif
 412       return
 413       end
 414 #endif
 415 c-----------------------------------------------------------------
 416       subroutine statout(itime)
 417       implicit real*8 (a-h,o-z)
 418       include 'DIMENSIONS'
 419       include 'COMMON.CONTROL'
 420       include 'COMMON.CHAIN'
 421       include 'COMMON.INTERACT'
 422       include 'COMMON.NAMES'
 423       include 'COMMON.IOUNITS'
 424       include 'COMMON.HEADER'
 425       include 'COMMON.SBRIDGE'
 426       include 'COMMON.DISTFIT'
 427       include 'COMMON.MD'
 428       include 'COMMON.REMD'
 429       include 'COMMON.SETUP'
 430       integer itime
 431       double precision energia(0:n_ene)
 432       double precision gyrate
 433       external gyrate
 434       common /gucio/ cm
 435       character*256 line1,line2
 436       character*4 format1,format2
 437       character*30 format
 438 #ifdef AIX
 439       if(itime.eq.0) then
 440        open(istat,file=statname,position="append")
 441       endif
 442 #else
 443 #ifdef PGI
 444       open(istat,file=statname,position="append")
 445 #else
 446       open(istat,file=statname,access="append")
 447 #endif
 448 #endif
 449        if (AFMlog.gt.0) then
 450        if (refstr) then
 451          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 452           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')
 453      &          itime,totT,EK,potE,totE,
 454      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
 455      &          potEcomp(23),me
 456           format1="a133"
 457          else
 458 C          print *,'A CHUJ',potEcomp(23)
 459           write (line1,'(i10,f15.2,7f12.3,i5,$)')
 460      &           itime,totT,EK,potE,totE,
 461      &           kinetic_T,t_bath,gyrate(),
 462      &           potEcomp(23),me
 463           format1="a114"
 464         endif
 465        else if (selfguide.gt.0) then
 466        distance=0.0
 467        do j=1,3
 468        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
 469        enddo
 470        distance=dsqrt(distance)
 471        if (refstr) then
 472          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 473           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2,
 474      &    f9.2,i5,$)')
 475      &          itime,totT,EK,potE,totE,
 476      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
 477      &          distance,potEcomp(23),me
 478           format1="a133"
 479 C          print *,"CHUJOWO"
 480          else
 481 C          print *,'A CHUJ',potEcomp(23)
 482           write (line1,'(i10,f15.2,8f12.3,i5,$)')
 483      &           itime,totT,EK,potE,totE,
 484      &           kinetic_T,t_bath,gyrate(),
 485      &           distance,potEcomp(23),me
 486           format1="a114"
 487         endif
 488        else
 489        if (refstr) then
 490          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 491           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')
 492      &          itime,totT,EK,potE,totE,
 493      &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 494           format1="a133"
 495         else
 496           write (line1,'(i10,f15.2,7f12.3,i5,$)')
 497      &           itime,totT,EK,potE,totE,
 498      &           amax,kinetic_T,t_bath,gyrate(),me
 499           format1="a114"
 500         endif
 501         endif
 502         if(usampl.and.totT.gt.eq_time) then
 503            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
 504      &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
 505      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 506            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
 507      &             +21*nfrag_back
 508         else
 509            format2="a001"
 510            line2=' '
 511         endif
 512         if (print_compon) then
 513           if(itime.eq.0) then
 514            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
 515      &                                                     ",20a12)"
 516            write (istat,format) "#","",
 517      &      (ename(print_order(i)),i=1,nprint_ene)
 518           endif
 519           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
 520      &                                                     ",20f12.3)"
 521           write (istat,format) line1,line2,
 522      &      (potEcomp(print_order(i)),i=1,nprint_ene)
 523         else
 524           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 525           write (istat,format) line1,line2
 526         endif
 527 #if defined(AIX)
 528         call flush(istat)
 529 #else
 530         close(istat)
 531 #endif
 532        return
 533       end
 534 c---------------------------------------------------------------
 535       double precision function gyrate()
 536       implicit real*8 (a-h,o-z)
 537       include 'DIMENSIONS'
 538       include 'COMMON.INTERACT'
 539       include 'COMMON.CHAIN'
 540       double precision cen(3),rg
 541
 542       do j=1,3
 543        cen(j)=0.0d0
 544       enddo
 545
 546       do i=nnt,nct
 547           do j=1,3
 548             cen(j)=cen(j)+c(j,i)
 549           enddo
 550       enddo
 551       do j=1,3
 552             cen(j)=cen(j)/dble(nct-nnt+1)
 553       enddo
 554       rg = 0.0d0
 555       do i = nnt, nct
 556         do j=1,3
 557          rg = rg + (c(j,i)-cen(j))**2
 558         enddo
 559       end do
 560       gyrate = sqrt(rg/dble(nct-nnt+1))
 561       return
 562       end
 563