parameters: (self, tol=DEFAULT_SYM_TOL)
return_statement: return frame
docstring: Determine symmetry reference frame. If noreorient is not set,
          this is the rotation matrix applied to the geometry in update_geometry.
  
          >>> print(H2OH2O.symmetry_frame())
          [[1.0, -0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -0.0, 1.0]]
docstring_summary: Determine symmetry reference frame. If noreorient is not set,
          this is the rotation matrix applied to the geometry in update_geometry.
docstring_tokens: Determine symmetry reference frame . If noreorient is not set this is the rotation matrix applied to the geometry in update_geometry .
function: def symmetry_frame(self, tol=DEFAULT_SYM_TOL):
          """Determine symmetry reference frame. If noreorient is not set,
          this is the rotation matrix applied to the geometry in update_geometry.
  
          >>> print(H2OH2O.symmetry_frame())
          [[1.0, -0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -0.0, 1.0]]
  
          """
          com = self.center_of_mass()
          if self.wholegeom is not None:
              current_geom = self.wholegeom
          else:
              current_geom = self.geometry(np_out=True)
          shifted_geom = current_geom - np.asarray(com)
          shifted_geom = shifted_geom.tolist()
  
          worldxaxis = [1.0, 0.0, 0.0]
          worldyaxis = [0.0, 1.0, 0.0]
          worldzaxis = [0.0, 0.0, 1.0]
  
          sigma = [0.0, 0.0, 0.0]
          sigmav = [0.0, 0.0, 0.0]
          c2axis = [0.0, 0.0, 0.0]
          c2axisperp = [0.0, 0.0, 0.0]
  
          linear, planar = self.is_linear_planar(tol)
          have_inversion = self.has_inversion(com, tol)
  
          # check for C2 axis
          have_c2axis = False
          if self.natom() < 2:
              have_c2axis = True
              c2axis = [0.0, 0.0, 1.0]
  
          elif linear:
              have_c2axis = True
              c2axis = sub(self.xyz(1), self.xyz(0))
              c2axis = normalize(c2axis)
  
          elif planar and have_inversion:
              # there is a c2 axis that won't be found using the usual
              #   algorithm. find two noncolinear atom-atom vectors (we know
              #   that linear == 0)
              BA = sub(self.xyz(1), self.xyz(0))
              BA = normalize(BA)
              for i in range(2, self.natom()):
                  CA = sub(self.xyz(i), self.xyz(0))
                  CA = normalize(CA)
                  BAxCA = cross(BA, CA)
                  if norm(BAxCA) > tol:
                      have_c2axis = True
                      BAxCA = normalize(BAxCA)
                      c2axis = copy.deepcopy(BAxCA)
                      break
  
          else:
              # loop through pairs of atoms to find c2 axis candidates
              for i in range(self.natom()):
                  A = shifted_geom[i]
                  AdotA = dot(A, A)
                  for j in range(i + 1):
                      # the atoms must be identical
                      if not self.atoms[i].is_equivalent_to(self.atoms[j]):
                          continue
                      B = shifted_geom[j]
                      # the atoms must be the same distance from the com
                      if math.fabs(AdotA - dot(B, B)) > tol:
                          continue
                      axis = add(A, B)
                      # atoms colinear with the com don't work
                      if norm(axis) < tol:
                          continue
                      axis = normalize(axis)
                      if self.is_axis(com, axis, 2, tol):
                          have_c2axis = True
                          c2axis = copy.deepcopy(axis)
                          break
                  else:
                      continue
                  break
  
          # symmframe found c2axis
          c2like = 'ZAxis'
          if have_c2axis:
              # try to make the sign of the axis correspond to one of the world axes
              c2like, c2axis = self.like_world_axis(c2axis, worldxaxis, worldyaxis, worldzaxis)
  
          # check for c2 axis perp to first c2 axis
          have_c2axisperp = False
          if have_c2axis:
              if self.natom() < 2:
                  have_c2axisperp = True
                  c2axisperp = [1.0, 0.0, 0.0]
  
              elif linear:
                  if have_inversion:
                      have_c2axisperp = True
                      c2axisperp = perp_unit(c2axis, [0.0, 0.0, 1.0])
  
              else:
                  # loop through paris of atoms to find c2 axis candidates
                  for i in range(self.natom()):
                      A = sub(self.xyz(i), com)
                      AdotA = dot(A, A)
                      for j in range(i):
                          # the atoms must be identical
                          if not self.atoms[i].is_equivalent_to(self.atoms[j]):
                              continue
                          B = sub(self.xyz(j), com)
                          # the atoms must be the same distance from the com
                          if math.fabs(AdotA - dot(B, B)) > tol:
                              continue
                          axis = add(A, B)
                          # atoms colinear with the com don't work
                          if norm(axis) < tol:
                              continue
                          axis = normalize(axis)
                          # if axis is not perp continue
                          if math.fabs(dot(axis, c2axis)) > tol:
                              continue
                          if self.is_axis(com, axis, 2, tol):
                              have_c2axisperp = True
                              c2axisperp = copy.deepcopy(axis)
                          break
                      else:
                          continue
                      break
  
          # symmframe found c2axisperp
          if have_c2axisperp:
              # try to make the sign of the axis correspond to one of the world axes
              c2perplike, c2axisperp = self.like_world_axis(c2axisperp, worldxaxis, worldyaxis, worldzaxis)
  
              # try to make c2axis the z axis
              if c2perplike == 'ZAxis':
                  tmpv = copy.deepcopy(c2axisperp)
                  c2axisperp = copy.deepcopy(c2axis)
                  c2axis = copy.deepcopy(tmpv)
                  c2perplike = c2like
                  c2like = 'ZAxis'
  
              if c2like != 'ZAxis':
                  if c2like == 'XAxis':
                      c2axis = cross(c2axis, c2axisperp)
                  else:
                      c2axis = cross(c2axisperp, c2axis)
                  c2like, c2axis = self.like_world_axis(c2axis, worldxaxis, worldyaxis, worldzaxis)
  
              # try to make c2axisperplike the x axis
              if c2perplike == 'YAxis':
                  c2axisperp = cross(c2axisperp, c2axis)
                  c2perplike, c2axisperp = self.like_world_axis(c2axisperp, worldxaxis, worldyaxis, worldzaxis)
  
          # Check for vertical plane
          have_sigmav = False
          if have_c2axis:
              if self.natom() < 2:
                  have_sigmav = True
                  sigmav = copy.deepcopy(c2axisperp)
  
              elif linear:
                  have_sigmav = True
                  if have_c2axisperp:
                      sigmav = copy.deepcopy(c2axisperp)
  
                  else:
                      sigmav = perp_unit(c2axis, [0.0, 0.0, 1.0])
              else:
                  # loop through pairs of atoms to find sigma v plane candidates
                  for i in range(self.natom()):
                      A = sub(self.xyz(i), com)
                      AdotA = dot(A, A)
                      # the second atom can equal i because i might be in the plane
                      for j in range(i + 1):
                          # the atoms must be identical
                          if not self.atoms[i].is_equivalent_to(self.atoms[j]):
                              continue
                          B = sub(self.xyz(j), com)
                          # the atoms must be the same distance from the com
                          if math.fabs(AdotA - dot(B, B)) > tol:
                              continue
                          inplane = add(B, A)
                          norm_inplane = norm(inplane)
                          if norm_inplane < tol:
                              continue
                          inplane = scale(inplane, 1.0 / norm_inplane)
                          perp = cross(c2axis, inplane)
                          norm_perp = norm(perp)
                          if norm_perp < tol:
                              continue
                          perp = scale(perp, 1.0 / norm_perp)
                          if self.is_plane(com, perp, tol):
                              have_sigmav = True
                              sigmav = copy.deepcopy(perp)
                              break
                      else:
                          continue
                      break
  
          # symmframe found sigmav
          if have_sigmav:
              # try to make the sign of the oop vec correspond to one of the world axes
              sigmavlike, sigmav = self.like_world_axis(sigmav, worldxaxis, worldyaxis, worldzaxis)
  
              # Choose sigmav to be the world x axis, if possible
              if c2like == 'ZAxis' and sigmavlike == 'YAxis':
                  sigmav = cross(sigmav, c2axis)
              elif c2like == 'YAxis' and sigmavlike == 'ZAxis':
                  sigmav = cross(c2axis, sigmav)
  
          # under certain conditions i need to know if there is any sigma plane
          have_sigma = False
          if not have_inversion and not have_c2axis:
              if planar:
                  # find two noncolinear atom-atom vectors
                  # we know that linear==0 since !have_c2axis
                  BA = sub(self.xyz(1), self.xyz(0))
                  BA = normalize(BA)
                  for i in range(2, self.natom()):
                      CA = sub(self.xyz(i), self.xyz(0))
                      CA = normalize(CA)
                      BAxCA = cross(BA, CA)
                      if norm(BAxCA) > tol:
                          have_sigma = True
                          BAxCA = normalize(BAxCA)
                          sigma = copy.deepcopy(BAxCA)
                          break
              else:
                  # loop through pairs of atoms to contruct trial planes
                  for i in range(self.natom()):
                      A = sub(self.xyz(i), com)
                      AdotA = dot(A, A)
                      for j in range(i):
                          # the atoms must be identical
                          if not self.atoms[i].is_equivalent_to(self.atoms[j]):
                              continue
                          B = sub(self.xyz(j), com)
                          BdotB = dot(B, B)
                          # the atoms must be the same distance from the com
                          if math.fabs(AdotA - BdotB) > tol:
                              continue
                          perp = sub(B, A)
                          norm_perp = norm(perp)
                          if norm_perp < tol:
                              continue
                          perp = scale(perp, 1.0 / norm_perp)
                          if self.is_plane(com, perp, tol):
                              have_sigma = True
                              sigma = copy.deepcopy(perp)
                              break
                      else:
                          continue
                      break
  
          # foundsigma
          if have_sigma:
              # try to make the sign of the oop vec correspond to one of the world axes
              xlikeness = math.fabs(dot(sigma, worldxaxis))
              ylikeness = math.fabs(dot(sigma, worldyaxis))
              zlikeness = math.fabs(dot(sigma, worldzaxis))
  
              if xlikeness > ylikeness and xlikeness > zlikeness:
                  if dot(sigma, worldxaxis) < 0:
                      sigma = scale(sigma, -1.0)
              elif ylikeness > zlikeness:
                  if dot(sigma, worldyaxis) < 0:
                      sigma = scale(sigma, -1.0)
              else:
                  if dot(sigma, worldzaxis) < 0:
                      sigma = scale(sigma, -1.0)
  
          # Find the three axes for the symmetry frame
          xaxis = copy.deepcopy(worldxaxis)
          zaxis = copy.deepcopy(worldzaxis)
          if have_c2axis:
              zaxis = copy.deepcopy(c2axis)
              if have_sigmav:
                  xaxis = copy.deepcopy(sigmav)
              elif have_c2axisperp:
                  xaxis = copy.deepcopy(c2axisperp)
              else:
                  # any axis orthogonal to the zaxis will do
                  xaxis = perp_unit(zaxis, zaxis)
          elif have_sigma:
              zaxis = copy.deepcopy(sigma)
              xaxis = perp_unit(zaxis, zaxis)
  
          # Clean up our z axis
          if math.fabs(zaxis[0]) < NOISY_ZERO:
              zaxis[0] = 0.0
          if math.fabs(zaxis[1]) < NOISY_ZERO:
              zaxis[1] = 0.0
          if math.fabs(zaxis[2]) < NOISY_ZERO:
              zaxis[2] = 0.0
  
          # Clean up our x axis
          if math.fabs(xaxis[0]) < NOISY_ZERO:
              xaxis[0] = 0.0
          if math.fabs(xaxis[1]) < NOISY_ZERO:
              xaxis[1] = 0.0
          if math.fabs(xaxis[2]) < NOISY_ZERO:
              xaxis[2] = 0.0
  
          # the y is then -x cross z
          yaxis = scale(cross(xaxis, zaxis), -1.0)
  
          #print("xaxis %20.14lf %20.14lf %20.14lf" % (xaxis[0], xaxis[1], xaxis[2]))
          #print("yaxis %20.14lf %20.14lf %20.14lf" % (yaxis[0], yaxis[1], yaxis[2]))
          #print("zaxis %20.14lf %20.14lf %20.14lf" % (zaxis[0], zaxis[1], zaxis[2]))
  
          frame = zero(3, 3)
          for i in range(3):
              frame[i][0] = xaxis[i]
              frame[i][1] = yaxis[i]
              frame[i][2] = zaxis[i]
          return frame
function_tokens: def symmetry_frame ( self , tol = DEFAULT_SYM_TOL ) : com = self . center_of_mass ( ) if self . wholegeom is not None : current_geom = self . wholegeom else : current_geom = self . geometry ( np_out = True ) shifted_geom = current_geom - np . asarray ( com ) shifted_geom = shifted_geom . tolist ( ) worldxaxis = [ 1.0 , 0.0 , 0.0 ] worldyaxis = [ 0.0 , 1.0 , 0.0 ] worldzaxis = [ 0.0 , 0.0 , 1.0 ] sigma = [ 0.0 , 0.0 , 0.0 ] sigmav = [ 0.0 , 0.0 , 0.0 ] c2axis = [ 0.0 , 0.0 , 0.0 ] c2axisperp = [ 0.0 , 0.0 , 0.0 ] linear , planar = self . is_linear_planar ( tol ) have_inversion = self . has_inversion ( com , tol ) # check for C2 axis have_c2axis = False if self . natom ( ) < 2 : have_c2axis = True c2axis = [ 0.0 , 0.0 , 1.0 ] elif linear : have_c2axis = True c2axis = sub ( self . xyz ( 1 ) , self . xyz ( 0 ) ) c2axis = normalize ( c2axis ) elif planar and have_inversion : # there is a c2 axis that won't be found using the usual #   algorithm. find two noncolinear atom-atom vectors (we know #   that linear == 0) BA = sub ( self . xyz ( 1 ) , self . xyz ( 0 ) ) BA = normalize ( BA ) for i in range ( 2 , self . natom ( ) ) : CA = sub ( self . xyz ( i ) , self . xyz ( 0 ) ) CA = normalize ( CA ) BAxCA = cross ( BA , CA ) if norm ( BAxCA ) > tol : have_c2axis = True BAxCA = normalize ( BAxCA ) c2axis = copy . deepcopy ( BAxCA ) break else : # loop through pairs of atoms to find c2 axis candidates for i in range ( self . natom ( ) ) : A = shifted_geom [ i ] AdotA = dot ( A , A ) for j in range ( i + 1 ) : # the atoms must be identical if not self . atoms [ i ] . is_equivalent_to ( self . atoms [ j ] ) : continue B = shifted_geom [ j ] # the atoms must be the same distance from the com if math . fabs ( AdotA - dot ( B , B ) ) > tol : continue axis = add ( A , B ) # atoms colinear with the com don't work if norm ( axis ) < tol : continue axis = normalize ( axis ) if self . is_axis ( com , axis , 2 , tol ) : have_c2axis = True c2axis = copy . deepcopy ( axis ) break else : continue break # symmframe found c2axis c2like = 'ZAxis' if have_c2axis : # try to make the sign of the axis correspond to one of the world axes c2like , c2axis = self . like_world_axis ( c2axis , worldxaxis , worldyaxis , worldzaxis ) # check for c2 axis perp to first c2 axis have_c2axisperp = False if have_c2axis : if self . natom ( ) < 2 : have_c2axisperp = True c2axisperp = [ 1.0 , 0.0 , 0.0 ] elif linear : if have_inversion : have_c2axisperp = True c2axisperp = perp_unit ( c2axis , [ 0.0 , 0.0 , 1.0 ] ) else : # loop through paris of atoms to find c2 axis candidates for i in range ( self . natom ( ) ) : A = sub ( self . xyz ( i ) , com ) AdotA = dot ( A , A ) for j in range ( i ) : # the atoms must be identical if not self . atoms [ i ] . is_equivalent_to ( self . atoms [ j ] ) : continue B = sub ( self . xyz ( j ) , com ) # the atoms must be the same distance from the com if math . fabs ( AdotA - dot ( B , B ) ) > tol : continue axis = add ( A , B ) # atoms colinear with the com don't work if norm ( axis ) < tol : continue axis = normalize ( axis ) # if axis is not perp continue if math . fabs ( dot ( axis , c2axis ) ) > tol : continue if self . is_axis ( com , axis , 2 , tol ) : have_c2axisperp = True c2axisperp = copy . deepcopy ( axis ) break else : continue break # symmframe found c2axisperp if have_c2axisperp : # try to make the sign of the axis correspond to one of the world axes c2perplike , c2axisperp = self . like_world_axis ( c2axisperp , worldxaxis , worldyaxis , worldzaxis ) # try to make c2axis the z axis if c2perplike == 'ZAxis' : tmpv = copy . deepcopy ( c2axisperp ) c2axisperp = copy . deepcopy ( c2axis ) c2axis = copy . deepcopy ( tmpv ) c2perplike = c2like c2like = 'ZAxis' if c2like != 'ZAxis' : if c2like == 'XAxis' : c2axis = cross ( c2axis , c2axisperp ) else : c2axis = cross ( c2axisperp , c2axis ) c2like , c2axis = self . like_world_axis ( c2axis , worldxaxis , worldyaxis , worldzaxis ) # try to make c2axisperplike the x axis if c2perplike == 'YAxis' : c2axisperp = cross ( c2axisperp , c2axis ) c2perplike , c2axisperp = self . like_world_axis ( c2axisperp , worldxaxis , worldyaxis , worldzaxis ) # Check for vertical plane have_sigmav = False if have_c2axis : if self . natom ( ) < 2 : have_sigmav = True sigmav = copy . deepcopy ( c2axisperp ) elif linear : have_sigmav = True if have_c2axisperp : sigmav = copy . deepcopy ( c2axisperp ) else : sigmav = perp_unit ( c2axis , [ 0.0 , 0.0 , 1.0 ] ) else : # loop through pairs of atoms to find sigma v plane candidates for i in range ( self . natom ( ) ) : A = sub ( self . xyz ( i ) , com ) AdotA = dot ( A , A ) # the second atom can equal i because i might be in the plane for j in range ( i + 1 ) : # the atoms must be identical if not self . atoms [ i ] . is_equivalent_to ( self . atoms [ j ] ) : continue B = sub ( self . xyz ( j ) , com ) # the atoms must be the same distance from the com if math . fabs ( AdotA - dot ( B , B ) ) > tol : continue inplane = add ( B , A ) norm_inplane = norm ( inplane ) if norm_inplane < tol : continue inplane = scale ( inplane , 1.0 / norm_inplane ) perp = cross ( c2axis , inplane ) norm_perp = norm ( perp ) if norm_perp < tol : continue perp = scale ( perp , 1.0 / norm_perp ) if self . is_plane ( com , perp , tol ) : have_sigmav = True sigmav = copy . deepcopy ( perp ) break else : continue break # symmframe found sigmav if have_sigmav : # try to make the sign of the oop vec correspond to one of the world axes sigmavlike , sigmav = self . like_world_axis ( sigmav , worldxaxis , worldyaxis , worldzaxis ) # Choose sigmav to be the world x axis, if possible if c2like == 'ZAxis' and sigmavlike == 'YAxis' : sigmav = cross ( sigmav , c2axis ) elif c2like == 'YAxis' and sigmavlike == 'ZAxis' : sigmav = cross ( c2axis , sigmav ) # under certain conditions i need to know if there is any sigma plane have_sigma = False if not have_inversion and not have_c2axis : if planar : # find two noncolinear atom-atom vectors # we know that linear==0 since !have_c2axis BA = sub ( self . xyz ( 1 ) , self . xyz ( 0 ) ) BA = normalize ( BA ) for i in range ( 2 , self . natom ( ) ) : CA = sub ( self . xyz ( i ) , self . xyz ( 0 ) ) CA = normalize ( CA ) BAxCA = cross ( BA , CA ) if norm ( BAxCA ) > tol : have_sigma = True BAxCA = normalize ( BAxCA ) sigma = copy . deepcopy ( BAxCA ) break else : # loop through pairs of atoms to contruct trial planes for i in range ( self . natom ( ) ) : A = sub ( self . xyz ( i ) , com ) AdotA = dot ( A , A ) for j in range ( i ) : # the atoms must be identical if not self . atoms [ i ] . is_equivalent_to ( self . atoms [ j ] ) : continue B = sub ( self . xyz ( j ) , com ) BdotB = dot ( B , B ) # the atoms must be the same distance from the com if math . fabs ( AdotA - BdotB ) > tol : continue perp = sub ( B , A ) norm_perp = norm ( perp ) if norm_perp < tol : continue perp = scale ( perp , 1.0 / norm_perp ) if self . is_plane ( com , perp , tol ) : have_sigma = True sigma = copy . deepcopy ( perp ) break else : continue break # foundsigma if have_sigma : # try to make the sign of the oop vec correspond to one of the world axes xlikeness = math . fabs ( dot ( sigma , worldxaxis ) ) ylikeness = math . fabs ( dot ( sigma , worldyaxis ) ) zlikeness = math . fabs ( dot ( sigma , worldzaxis ) ) if xlikeness > ylikeness and xlikeness > zlikeness : if dot ( sigma , worldxaxis ) < 0 : sigma = scale ( sigma , - 1.0 ) elif ylikeness > zlikeness : if dot ( sigma , worldyaxis ) < 0 : sigma = scale ( sigma , - 1.0 ) else : if dot ( sigma , worldzaxis ) < 0 : sigma = scale ( sigma , - 1.0 ) # Find the three axes for the symmetry frame xaxis = copy . deepcopy ( worldxaxis ) zaxis = copy . deepcopy ( worldzaxis ) if have_c2axis : zaxis = copy . deepcopy ( c2axis ) if have_sigmav : xaxis = copy . deepcopy ( sigmav ) elif have_c2axisperp : xaxis = copy . deepcopy ( c2axisperp ) else : # any axis orthogonal to the zaxis will do xaxis = perp_unit ( zaxis , zaxis ) elif have_sigma : zaxis = copy . deepcopy ( sigma ) xaxis = perp_unit ( zaxis , zaxis ) # Clean up our z axis if math . fabs ( zaxis [ 0 ] ) < NOISY_ZERO : zaxis [ 0 ] = 0.0 if math . fabs ( zaxis [ 1 ] ) < NOISY_ZERO : zaxis [ 1 ] = 0.0 if math . fabs ( zaxis [ 2 ] ) < NOISY_ZERO : zaxis [ 2 ] = 0.0 # Clean up our x axis if math . fabs ( xaxis [ 0 ] ) < NOISY_ZERO : xaxis [ 0 ] = 0.0 if math . fabs ( xaxis [ 1 ] ) < NOISY_ZERO : xaxis [ 1 ] = 0.0 if math . fabs ( xaxis [ 2 ] ) < NOISY_ZERO : xaxis [ 2 ] = 0.0 # the y is then -x cross z yaxis = scale ( cross ( xaxis , zaxis ) , - 1.0 ) #print("xaxis %20.14lf %20.14lf %20.14lf" % (xaxis[0], xaxis[1], xaxis[2])) #print("yaxis %20.14lf %20.14lf %20.14lf" % (yaxis[0], yaxis[1], yaxis[2])) #print("zaxis %20.14lf %20.14lf %20.14lf" % (zaxis[0], zaxis[1], zaxis[2])) frame = zero ( 3 , 3 ) for i in range ( 3 ) : frame [ i ] [ 0 ] = xaxis [ i ] frame [ i ] [ 1 ] = yaxis [ i ] frame [ i ] [ 2 ] = zaxis [ i ] return frame