c                        PROGRAM CENTROID           E. Meng  1/94
c
c     calculates centroid and radius of gyration (both NON-MASS-WEIGHTED)
c     given an input PDB file
c
      integer maxats
      parameter (maxats=10000)
      real ca(maxats,3), c(3), com(3), invna
      character*4 aresa(maxats), ares, dum
      character*4 atypa(maxats), atyp
      character*4 inuma(maxats), inum
      character*4 irnma(maxats), irnm
      character*80 aname
      integer i, j, na
      real dsq, sumdsq, rogsq
c
      write (6, *) 'input PDB filename?'
      read (5, 1000) aname
 1000 format (A80)
c
c  --read in the PDB file
c
      open (unit=1, file=aname, status='old')
      i=0
   10 read (1, 1002, end=50) dum, inum, atyp, ares, irnm,
     &(c(j), j=1, 3)
 1002 format (A4, 3x, 3(A4, 1x), A4, 4x, 3F8.3)
      if ((dum .ne. 'ATOM').and.(dum .ne. 'HETA')) go to 10
      i=i + 1
      if (i .gt. maxats) then
        write (6, *) ' atom limit exceeded--program stops'
        stop
      endif
      atypa(i)=atyp
      irnma(i)=irnm
      aresa(i)=ares
      inuma(i)=inum
      do 20 j=1, 3
        ca(i,j)=c(j)
   20 continue
      go to 10
   50 write (6, '(I5, A19)') i, ' atoms in structure'
      na=i
      close (1)
c
c  --calculate the non-mass-weighted centroid
c
      invna=1.0/float(na)
      do 100 j=1, 3
        com(j)=0.0
  100 continue
      do 120 i=1, na
        do 110 j=1, 3
          com(j)=com(j) + ca(i,j)
  110   continue
  120 continue
      do 130 j=1, 3
        com(j)=com(j)*invna
  130 continue
      write (6,'(A, 3F8.3)') '  centroid: ',com(1),com(2),com(3)
c
c  --calculate the non-mass-weighted radius of gyration
c
      sumdsq = 0.0
      do 200 i=1, na
        dsq = 0.0
        do 140 j=1, 3
          dsq = dsq + (ca(i,j)-com(j))**2
  140   continue
        sumdsq = sumdsq + dsq
  200 continue
      rogsq = sumdsq*invna
      write (6,*) ' mean square distance from centroid = ', rogsq
      write (6,*) ' radius of gyration = ', sqrt(rogsq)
      end