subroutine cmshft(natom,rx,ry,rz,rmass,cmx,cmy,cmz,ishift)
  !c
  !c      shifts all x, y, and z coordinates passed in rx,ry,rz to keep
  !c      the center of mass at the origin if ishift is not zero. otherwise
  !c      cmx, cmy, and cmz return the center of mass coordinates and no
  !c      shift is performed. coordinate and mass units arbitrary.
  !c

  use VarTypes
  use Io_module

  IMPLICIT NONE

  
! Number of atoms
  INTEGER(KINT), INTENT(IN) :: NAtom
! Coordinates and masses 
  REAL(KREAL), DIMENSION(NAtom) :: rx,ry,rz,rmass
! Center of mass coordinates
  real(KREAL) :: cmx,cmy,cmz
! Flag if user wants to shift center of mass to origin
  INTEGER(KINT) :: IShift

  REAL(KREAL), parameter :: zero=0.d0,one=1.d0
  INTEGER(KINT) :: j
  REAL(KREAL) :: TmassI, tmass

  if(natom.le.0) then
     if(IOOUT.gt.0) write(IOOUT,'(/2a/)')  &
          '  *** error, illegal array dimension in cmshft,', &
          ' return without action ***'
     return
  end if
  !c
  !c      calculate center of mass coordinates
  !c
  tmass=rmass(1)
  do j=2,natom
     tmass=tmass+rmass(j)
  END DO
  !c
  !c     DO i=1,natom
  !c       write(IOOUT,'(4F15.3)') rmass(i),rx(i),ry(i),rz(i)
  !c     END DO

  cmx=zero
  cmy=zero
  cmz=zero
  do j=1,natom
     cmx=rmass(j)*rx(j)+cmx
     cmy=rmass(j)*ry(j)+cmy
     cmz=rmass(j)*rz(j)+cmz
  END DO

  tmassi=one/tmass
  cmx=cmx*tmassi
  cmy=cmy*tmassi
  cmz=cmz*tmassi
  !c
  !c      optional shift of all cartesian coordinates
  !c
  if(ishift.ne.0) then
     do j=1,natom
        rx(j)=rx(j)-cmx
        ry(j)=ry(j)-cmy
        rz(j)=rz(j)-cmz
     END DO
  end if

  return
end subroutine cmshft
