code: def get_gyration_tensor(elements, coordinates):
      
      # First calculate COM for correction.
      com = centre_of_mass(elements, coordinates)
      # Correct the coordinates for the COM.
      coordinates = coordinates - com
      # Calculate diagonal and then other values of the matrix.
      diag = np.sum(coordinates**2, axis=0)
      xy = np.sum(coordinates[:, 0] * coordinates[:, 1])
      xz = np.sum(coordinates[:, 0] * coordinates[:, 2])
      yz = np.sum(coordinates[:, 1] * coordinates[:, 2])
      S = np.array([[diag[0], xy, xz], [xy, diag[1], yz],
                    [xz, yz, diag[2]]]) / coordinates.shape[0]
      return (S)
docs: Return the gyration tensor of a molecule.
  
      The gyration tensor should be invariant to the molecule's position.
      The known formulas for the gyration tensor have the correction for the
      centre of mass of the molecule, therefore, the coordinates are first
      corrected for the centre of mass and essentially shifted to the origin.
  
      Parameters
      ----------
      elements : numpy.ndarray
          The array containing the molecule's elemental data.
  
      coordinates :  numpy.ndarray
          The array containing the Cartesian coordinates of the molecule.
  
      Returns
      -------
      numpy.ndarray
          The gyration tensor of a molecule invariant to the molecule's position.