parameters: (self, tol=FULL_PG_TOL)
return_statement: return
docstring: Determine and set FULL point group. self.PYfull_pg_n is highest
          order n in Cn. 0 for atoms or infinity.
docstring_summary: Determine and set FULL point group. self.PYfull_pg_n is highest
          order n in Cn. 0 for atoms or infinity.
docstring_tokens: Determine and set FULL point group . self . PYfull_pg_n is highest order n in Cn . 0 for atoms or infinity .
function: def set_full_point_group(self, tol=FULL_PG_TOL):
          """Determine and set FULL point group. self.PYfull_pg_n is highest
          order n in Cn. 0 for atoms or infinity.
  
          """
          verbose = 1  # TODO
          # Get cartesian geometry and put COM at origin
          geom = self.geometry()
          com = self.center_of_mass()
          for at in range(self.natom()):
              geom[at][0] += -com[0]
              geom[at][1] += -com[1]
              geom[at][2] += -com[2]
  
          # Get rotor type
          rotor = self.rotor_type(tol)
          if verbose > 2:
              print("""  Rotor type                       : %s""" % (rotor))
  
          # Get the D2h point group from Jet and Ed's code: c1 ci c2 cs d2 c2v c2h d2h
          #   and ignore the user-specified subgroup in this case.
          pg = self.find_highest_point_group(tol)
          d2h_subgroup = pg.symbol()
          if verbose > 2:
              print("""  D2h_subgroup                     : %s""" % (self.point_group().symbol()))
  
          # Check inversion
          v3_zero = [0.0, 0.0, 0.0]
          op_i = self.has_inversion(v3_zero, tol)
          if verbose > 2:
              print("""  Inversion symmetry               : %s""" % ('yes' if op_i else 'no'))
  
          x_axis = [1, 0, 0]
          y_axis = [0, 1, 0]
          z_axis = [0, 0, 1]
          rot_axis = [0.0, 0.0, 0.0]
  
          if rotor == 'RT_ATOM':  # atoms
              self.full_pg = 'ATOM'
              self.PYfull_pg_n = 0
  
          elif rotor == 'RT_LINEAR':  # linear molecules
              self.full_pg = 'D_inf_h' if op_i else 'C_inf_v'
              self.PYfull_pg_n = 0
  
          elif rotor == 'RT_SPHERICAL_TOP':  # spherical tops
              if not op_i:  # The only spherical top without inversion is Td.
                  self.full_pg = 'Td'
                  self.PYfull_pg_n = 3
              else:  # Oh or Ih ?
                  # Oh has a S4 and should be oriented properly already.
                  test_mat = matrix_3d_rotation(geom, z_axis, math.pi / 2.0, True)
                  op_symm = equal_but_for_row_order(geom, test_mat, tol)
                  if verbose > 2:
                      print("""  S4z                              : %s""" % ('yes' if op_symm else 'no'))
  
                  if op_symm:
                      self.full_pg = 'Oh'
                      self.PYfull_pg_n = 4
                  else:
                      self.full_pg = 'Ih'
                      self.PYfull_pg_n = 5
  
          elif rotor == 'RT_ASYMMETRIC_TOP':  # asymmetric tops cannot exceed D2h, right?
  
              if d2h_subgroup == 'c1':
                  self.full_pg = 'C1'
                  self.PYfull_pg_n = 1
  
              elif d2h_subgroup == 'ci':
                  self.full_pg = 'Ci'
                  self.PYfull_pg_n = 1
  
              elif d2h_subgroup == 'c2':
                  self.full_pg = 'Cn'
                  self.PYfull_pg_n = 2
  
              elif d2h_subgroup == 'cs':
                  self.full_pg = 'Cs'
                  self.PYfull_pg_n = 1
  
              elif d2h_subgroup == 'd2':
                  self.full_pg = 'Dn'
                  self.PYfull_pg_n = 2
  
              elif d2h_subgroup == 'c2v':
                  self.full_pg = 'Cnv'
                  self.PYfull_pg_n = 2
  
              elif d2h_subgroup == 'c2h':
                  self.full_pg = 'Cnh'
                  self.PYfull_pg_n = 2
  
              elif d2h_subgroup == 'd2h':
                  self.full_pg = 'Dnh'
                  self.PYfull_pg_n = 2
  
              else:
                  print("""  Warning: Cannot determine point group.""")
  
          elif rotor in ['RT_SYMMETRIC_TOP', 'RT_PROLATE_SYMMETRIC_TOP', 'RT_OBLATE_SYMMETRIC_TOP']:
  
              # Find principal axis that is unique and make it z-axis.
              It = self.inertia_tensor()
              I_evals, I_evecs = diagonalize3x3symmat(It)
              ev_list = list(zip(I_evals, transpose(I_evecs)))  # eigenvectors are cols of I_evecs
              ev_list.sort(key=lambda tup: tup[0], reverse=False)
              I_evals, I_evecs = zip(*ev_list)  # sorted eigenvectors are now rows of I_evecs
              if verbose > 2:
                  print("""  I_evals: %15.10lf %15.10lf %15.10lf""" % (I_evals[0], I_evals[1], I_evals[2]))
  
              unique_axis = 1
              if abs(I_evals[0] - I_evals[1]) < tol:
                  unique_axis = 2
              elif abs(I_evals[1] - I_evals[2]) < tol:
                  unique_axis = 0
  
              # Compute angle between unique axis and the z-axis
              old_axis = I_evecs[unique_axis]
  
              ddot = dot(z_axis, old_axis)
              if abs(ddot - 1) < 1.0e-10:
                  phi = 0.0
              elif abs(ddot + 1) < 1.0e-10:
                  phi = math.pi
              else:
                  phi = math.acos(ddot)
  
              # Rotate geometry to put unique axis on the z-axis, if it isn't already.
              if abs(phi) > 1.0e-14:
                  rot_axis = cross(z_axis, old_axis)  # right order?
                  test_mat = matrix_3d_rotation(geom, rot_axis, phi, False)
                  if verbose > 2:
                      print("""  Rotating by %lf to get principal axis on z-axis ...""" % (phi))
                  geom = [row[:] for row in test_mat]
  
              if verbose > 2:
                  print("""  Geometry to analyze - principal axis on z-axis:""")
                  for at in range(self.natom()):
                      print("""%20.15lf %20.15lf %20.15lf""" % (geom[at][0], geom[at][1], geom[at][2]))
                  print('\n')
  
              # Determine order Cn and Sn of principal axis.
              Cn_z = matrix_3d_rotation_Cn(geom, z_axis, False, tol)
              if verbose > 2:
                  print("""  Highest rotation axis (Cn_z)     : %d""" % (Cn_z))
  
              Sn_z = matrix_3d_rotation_Cn(geom, z_axis, True, tol)
              if verbose > 2:
                  print("""  Rotation axis (Sn_z)             : %d""" % (Sn_z))
  
              # Check for sigma_h (xy plane).
              op_sigma_h = False
              for at in range(self.natom()):
                  if abs(geom[at][2]) < tol:
                      continue  # atom is in xy plane
                  else:
                      test_atom = [geom[at][0], geom[at][1], -1 * geom[at][2]]
                      if not atom_present_in_geom(geom, test_atom, tol):
                          break
              else:
                  op_sigma_h = True
              if verbose > 2:
                  print("""  sigma_h                          : %s""" % ('yes' if op_sigma_h else 'no'))
  
              # Rotate one off-axis atom to the yz plane and check for sigma_v's.
              for at in range(self.natom()):
                  dist_from_z = math.sqrt(geom[at][0] * geom[at][0] + geom[at][1] * geom[at][1])
                  if abs(dist_from_z) > tol:
                      pivot_atom_i = at
                      break
  
              if pivot_atom_i == self.natom():  # needs to be in else clause?
                  raise ValidationError("Molecule::set_full_point_group: Not a linear molecule but could not find off-axis atom.")
  
              # Rotate around z-axis to put pivot atom in the yz plane
              xy_point = normalize([geom[pivot_atom_i][0], geom[pivot_atom_i][1], 0])
              ddot = dot(y_axis, xy_point)
              phi = None
              if abs(ddot - 1) < 1.0e-10:
                  phi = 0.0
              elif abs(ddot + 1) < 1.0e-10:
                  phi = math.pi
              else:
                  phi = math.acos(ddot)
  
              is_D = False
              if abs(phi) > 1.0e-14:
                  test_mat = matrix_3d_rotation(geom, z_axis, phi, False)
                  if verbose > 2:
                      print("""  Rotating by %8.3e to get atom %d in yz-plane ...""" % (phi, pivot_atom_i + 1))
                  geom = [row[:] for row in test_mat]
  
              # Check for sigma_v (yz plane).
              op_sigma_v = False
              for at in range(self.natom()):
                  if abs(geom[at][0]) < tol:
                      continue  # atom is in yz plane
                  else:
                      test_atom = [-1 * geom[at][0], geom[at][1], geom[at][2]]
                      if not atom_present_in_geom(geom, test_atom, tol):
                          break
              else:
                  op_sigma_v = True
              if verbose > 2:
                  print("""  sigma_v                          : %s""" % ('yes' if op_sigma_v else 'no'))
  
                  print("""  geom to analyze - one atom in yz plane:""")
                  for at in range(self.natom()):
                      print("""%20.15lf %20.15lf %20.15lf""" % (geom[at][0], geom[at][1], geom[at][2]))
                  print('\n')
  
              # Check for perpendicular C2's.
              # Loop through pairs of atoms to find c2 axis candidates.
              for i in range(self.natom()):
                  A = [geom[i][0], geom[i][1], geom[i][2]]
                  AdotA = dot(A, A)
                  for j in range(i):
                      if self.Z(at) != self.Z(j):
                          continue  # ensure same atomic number
  
                      B = [geom[j][0], geom[j][1], geom[j][2]]  # ensure same distance from com
                      if abs(AdotA - dot(B, B)) > 1.0e-6:
                          continue  # loose check
  
                      # Use sum of atom vectors as axis if not 0.
                      axis = add(A, B)
                      if norm(axis) < 1.0e-12:
                          continue
                      axis = normalize(axis)
  
                      # Check if axis is perpendicular to z-axis.
                      if abs(dot(axis, z_axis)) > 1.0e-6:
                          continue
  
                      # Do the thorough check for C2.
                      if matrix_3d_rotation_Cn(geom, axis, False, tol, 2) == 2:
                          is_D = True
              if verbose > 2:
                  print("""  perp. C2's                       : %s""" % ('yes' if is_D else 'no'))
  
              # Now assign point groups!  Sn first.
              if Sn_z == 2 * Cn_z and not is_D:
                  self.full_pg = 'Sn'
                  self.PYfull_pg_n = Sn_z
                  return
  
              if is_D:  # has perpendicular C2's
                  if op_sigma_h and op_sigma_v:  # Dnh : Cn, nC2, sigma_h, nSigma_v
                      self.full_pg = 'Dnh'
                      self.PYfull_pg_n = Cn_z
  
                  elif Sn_z == 2 * Cn_z:  # Dnd : Cn, nC2, S2n axis coincident with Cn
                      self.full_pg = 'Dnd'
                      self.PYfull_pg_n = Cn_z
  
                  else:  # Dn : Cn, nC2
                      self.full_pg = 'Dn'
                      self.PYfull_pg_n = Cn_z
  
              else:  # lacks perpendicular C2's
                  if op_sigma_h and Sn_z == Cn_z:  # Cnh : Cn, sigma_h, Sn coincident with Cn
                      self.full_pg = 'Cnh'
                      self.PYfull_pg_n = Cn_z
  
                  elif op_sigma_v:  # Cnv : Cn, nCv
                      self.full_pg = 'Cnv'
                      self.PYfull_pg_n = Cn_z
  
                  else:  # Cn  : Cn
                      self.full_pg = 'Cn'
                      self.PYfull_pg_n = Cn_z
          return
function_tokens: def set_full_point_group ( self , tol = FULL_PG_TOL ) : verbose = 1 # TODO # Get cartesian geometry and put COM at origin geom = self . geometry ( ) com = self . center_of_mass ( ) for at in range ( self . natom ( ) ) : geom [ at ] [ 0 ] += - com [ 0 ] geom [ at ] [ 1 ] += - com [ 1 ] geom [ at ] [ 2 ] += - com [ 2 ] # Get rotor type rotor = self . rotor_type ( tol ) if verbose > 2 : print ( """  Rotor type                       : %s""" % ( rotor ) ) # Get the D2h point group from Jet and Ed's code: c1 ci c2 cs d2 c2v c2h d2h #   and ignore the user-specified subgroup in this case. pg = self . find_highest_point_group ( tol ) d2h_subgroup = pg . symbol ( ) if verbose > 2 : print ( """  D2h_subgroup                     : %s""" % ( self . point_group ( ) . symbol ( ) ) ) # Check inversion v3_zero = [ 0.0 , 0.0 , 0.0 ] op_i = self . has_inversion ( v3_zero , tol ) if verbose > 2 : print ( """  Inversion symmetry               : %s""" % ( 'yes' if op_i else 'no' ) ) x_axis = [ 1 , 0 , 0 ] y_axis = [ 0 , 1 , 0 ] z_axis = [ 0 , 0 , 1 ] rot_axis = [ 0.0 , 0.0 , 0.0 ] if rotor == 'RT_ATOM' : # atoms self . full_pg = 'ATOM' self . PYfull_pg_n = 0 elif rotor == 'RT_LINEAR' : # linear molecules self . full_pg = 'D_inf_h' if op_i else 'C_inf_v' self . PYfull_pg_n = 0 elif rotor == 'RT_SPHERICAL_TOP' : # spherical tops if not op_i : # The only spherical top without inversion is Td. self . full_pg = 'Td' self . PYfull_pg_n = 3 else : # Oh or Ih ? # Oh has a S4 and should be oriented properly already. test_mat = matrix_3d_rotation ( geom , z_axis , math . pi / 2.0 , True ) op_symm = equal_but_for_row_order ( geom , test_mat , tol ) if verbose > 2 : print ( """  S4z                              : %s""" % ( 'yes' if op_symm else 'no' ) ) if op_symm : self . full_pg = 'Oh' self . PYfull_pg_n = 4 else : self . full_pg = 'Ih' self . PYfull_pg_n = 5 elif rotor == 'RT_ASYMMETRIC_TOP' : # asymmetric tops cannot exceed D2h, right? if d2h_subgroup == 'c1' : self . full_pg = 'C1' self . PYfull_pg_n = 1 elif d2h_subgroup == 'ci' : self . full_pg = 'Ci' self . PYfull_pg_n = 1 elif d2h_subgroup == 'c2' : self . full_pg = 'Cn' self . PYfull_pg_n = 2 elif d2h_subgroup == 'cs' : self . full_pg = 'Cs' self . PYfull_pg_n = 1 elif d2h_subgroup == 'd2' : self . full_pg = 'Dn' self . PYfull_pg_n = 2 elif d2h_subgroup == 'c2v' : self . full_pg = 'Cnv' self . PYfull_pg_n = 2 elif d2h_subgroup == 'c2h' : self . full_pg = 'Cnh' self . PYfull_pg_n = 2 elif d2h_subgroup == 'd2h' : self . full_pg = 'Dnh' self . PYfull_pg_n = 2 else : print ( """  Warning: Cannot determine point group.""" ) elif rotor in [ 'RT_SYMMETRIC_TOP' , 'RT_PROLATE_SYMMETRIC_TOP' , 'RT_OBLATE_SYMMETRIC_TOP' ] : # Find principal axis that is unique and make it z-axis. It = self . inertia_tensor ( ) I_evals , I_evecs = diagonalize3x3symmat ( It ) ev_list = list ( zip ( I_evals , transpose ( I_evecs ) ) ) # eigenvectors are cols of I_evecs ev_list . sort ( key = lambda tup : tup [ 0 ] , reverse = False ) I_evals , I_evecs = zip ( * ev_list ) # sorted eigenvectors are now rows of I_evecs if verbose > 2 : print ( """  I_evals: %15.10lf %15.10lf %15.10lf""" % ( I_evals [ 0 ] , I_evals [ 1 ] , I_evals [ 2 ] ) ) unique_axis = 1 if abs ( I_evals [ 0 ] - I_evals [ 1 ] ) < tol : unique_axis = 2 elif abs ( I_evals [ 1 ] - I_evals [ 2 ] ) < tol : unique_axis = 0 # Compute angle between unique axis and the z-axis old_axis = I_evecs [ unique_axis ] ddot = dot ( z_axis , old_axis ) if abs ( ddot - 1 ) < 1.0e-10 : phi = 0.0 elif abs ( ddot + 1 ) < 1.0e-10 : phi = math . pi else : phi = math . acos ( ddot ) # Rotate geometry to put unique axis on the z-axis, if it isn't already. if abs ( phi ) > 1.0e-14 : rot_axis = cross ( z_axis , old_axis ) # right order? test_mat = matrix_3d_rotation ( geom , rot_axis , phi , False ) if verbose > 2 : print ( """  Rotating by %lf to get principal axis on z-axis ...""" % ( phi ) ) geom = [ row [ : ] for row in test_mat ] if verbose > 2 : print ( """  Geometry to analyze - principal axis on z-axis:""" ) for at in range ( self . natom ( ) ) : print ( """%20.15lf %20.15lf %20.15lf""" % ( geom [ at ] [ 0 ] , geom [ at ] [ 1 ] , geom [ at ] [ 2 ] ) ) print ( '\n' ) # Determine order Cn and Sn of principal axis. Cn_z = matrix_3d_rotation_Cn ( geom , z_axis , False , tol ) if verbose > 2 : print ( """  Highest rotation axis (Cn_z)     : %d""" % ( Cn_z ) ) Sn_z = matrix_3d_rotation_Cn ( geom , z_axis , True , tol ) if verbose > 2 : print ( """  Rotation axis (Sn_z)             : %d""" % ( Sn_z ) ) # Check for sigma_h (xy plane). op_sigma_h = False for at in range ( self . natom ( ) ) : if abs ( geom [ at ] [ 2 ] ) < tol : continue # atom is in xy plane else : test_atom = [ geom [ at ] [ 0 ] , geom [ at ] [ 1 ] , - 1 * geom [ at ] [ 2 ] ] if not atom_present_in_geom ( geom , test_atom , tol ) : break else : op_sigma_h = True if verbose > 2 : print ( """  sigma_h                          : %s""" % ( 'yes' if op_sigma_h else 'no' ) ) # Rotate one off-axis atom to the yz plane and check for sigma_v's. for at in range ( self . natom ( ) ) : dist_from_z = math . sqrt ( geom [ at ] [ 0 ] * geom [ at ] [ 0 ] + geom [ at ] [ 1 ] * geom [ at ] [ 1 ] ) if abs ( dist_from_z ) > tol : pivot_atom_i = at break if pivot_atom_i == self . natom ( ) : # needs to be in else clause? raise ValidationError ( "Molecule::set_full_point_group: Not a linear molecule but could not find off-axis atom." ) # Rotate around z-axis to put pivot atom in the yz plane xy_point = normalize ( [ geom [ pivot_atom_i ] [ 0 ] , geom [ pivot_atom_i ] [ 1 ] , 0 ] ) ddot = dot ( y_axis , xy_point ) phi = None if abs ( ddot - 1 ) < 1.0e-10 : phi = 0.0 elif abs ( ddot + 1 ) < 1.0e-10 : phi = math . pi else : phi = math . acos ( ddot ) is_D = False if abs ( phi ) > 1.0e-14 : test_mat = matrix_3d_rotation ( geom , z_axis , phi , False ) if verbose > 2 : print ( """  Rotating by %8.3e to get atom %d in yz-plane ...""" % ( phi , pivot_atom_i + 1 ) ) geom = [ row [ : ] for row in test_mat ] # Check for sigma_v (yz plane). op_sigma_v = False for at in range ( self . natom ( ) ) : if abs ( geom [ at ] [ 0 ] ) < tol : continue # atom is in yz plane else : test_atom = [ - 1 * geom [ at ] [ 0 ] , geom [ at ] [ 1 ] , geom [ at ] [ 2 ] ] if not atom_present_in_geom ( geom , test_atom , tol ) : break else : op_sigma_v = True if verbose > 2 : print ( """  sigma_v                          : %s""" % ( 'yes' if op_sigma_v else 'no' ) ) print ( """  geom to analyze - one atom in yz plane:""" ) for at in range ( self . natom ( ) ) : print ( """%20.15lf %20.15lf %20.15lf""" % ( geom [ at ] [ 0 ] , geom [ at ] [ 1 ] , geom [ at ] [ 2 ] ) ) print ( '\n' ) # Check for perpendicular C2's. # Loop through pairs of atoms to find c2 axis candidates. for i in range ( self . natom ( ) ) : A = [ geom [ i ] [ 0 ] , geom [ i ] [ 1 ] , geom [ i ] [ 2 ] ] AdotA = dot ( A , A ) for j in range ( i ) : if self . Z ( at ) != self . Z ( j ) : continue # ensure same atomic number B = [ geom [ j ] [ 0 ] , geom [ j ] [ 1 ] , geom [ j ] [ 2 ] ] # ensure same distance from com if abs ( AdotA - dot ( B , B ) ) > 1.0e-6 : continue # loose check # Use sum of atom vectors as axis if not 0. axis = add ( A , B ) if norm ( axis ) < 1.0e-12 : continue axis = normalize ( axis ) # Check if axis is perpendicular to z-axis. if abs ( dot ( axis , z_axis ) ) > 1.0e-6 : continue # Do the thorough check for C2. if matrix_3d_rotation_Cn ( geom , axis , False , tol , 2 ) == 2 : is_D = True if verbose > 2 : print ( """  perp. C2's                       : %s""" % ( 'yes' if is_D else 'no' ) ) # Now assign point groups!  Sn first. if Sn_z == 2 * Cn_z and not is_D : self . full_pg = 'Sn' self . PYfull_pg_n = Sn_z return if is_D : # has perpendicular C2's if op_sigma_h and op_sigma_v : # Dnh : Cn, nC2, sigma_h, nSigma_v self . full_pg = 'Dnh' self . PYfull_pg_n = Cn_z elif Sn_z == 2 * Cn_z : # Dnd : Cn, nC2, S2n axis coincident with Cn self . full_pg = 'Dnd' self . PYfull_pg_n = Cn_z else : # Dn : Cn, nC2 self . full_pg = 'Dn' self . PYfull_pg_n = Cn_z else : # lacks perpendicular C2's if op_sigma_h and Sn_z == Cn_z : # Cnh : Cn, sigma_h, Sn coincident with Cn self . full_pg = 'Cnh' self . PYfull_pg_n = Cn_z elif op_sigma_v : # Cnv : Cn, nCv self . full_pg = 'Cnv' self . PYfull_pg_n = Cn_z else : # Cn  : Cn self . full_pg = 'Cn' self . PYfull_pg_n = Cn_z return