parameters: (findifrec, freq_irrep_only)
return_statement: return _process_hessian(H_pi, B_pi, massweighter, data["print_lvl"])
docstring: Compute the Hessian by finite difference of energies.
  
      Parameters
      ----------
      findifrec : dict
          Dictionary of finite difference data, specified in _geom_generator docstring.
      freq_irrep_only : int
          The 0-indexed Cotton ordered irrep to get frequencies for. Choose -1 for all irreps.
  
      Returns
      -------
      hessian : ndarray
          (3 * nat, 3 * nat) Cartesian Hessian [Eh/a0^2].
docstring_summary: Compute the Hessian by finite difference of energies.
docstring_tokens: Compute the Hessian by finite difference of energies .
function: def assemble_hessian_from_energies(findifrec, freq_irrep_only):
      """Compute the Hessian by finite difference of energies.
  
      Parameters
      ----------
      findifrec : dict
          Dictionary of finite difference data, specified in _geom_generator docstring.
      freq_irrep_only : int
          The 0-indexed Cotton ordered irrep to get frequencies for. Choose -1 for all irreps.
  
      Returns
      -------
      hessian : ndarray
          (3 * nat, 3 * nat) Cartesian Hessian [Eh/a0^2].
      """
  
      # This *must* be a Psi molecule at present - CdSalcList generation panics otherwise
      mol = core.Molecule.from_schema(findifrec["molecule"], nonphysical=True, verbose=0)
  
      displacements = findifrec["displacements"]
      ref_energy = findifrec["reference"]["energy"]
  
      def init_string(data):
          max_label_len = str(max([len(label) for label in displacements]))
          out_str = ""
          for label, disp_data in displacements.items():
              out_str += ("    {:" + max_label_len + "s} : {:20.10f}\n").format(label, disp_data["energy"])
          return ("  Computing second-derivative from energies using projected, \n"
                  "  symmetry-adapted, cartesian coordinates.\n\n"
                  "  {:d} energies passed in, including the reference geometry.\n"
                  "    Using {:d}-point formula.\n"
                  "    Energy without displacement: {:15.10f}\n"
                  "    Check energies below for precision!\n{}".format(
                      len(displacements) + 1, findifrec["stencil_size"], ref_energy, out_str))
  
      data = _initialize_findif(mol, freq_irrep_only, "2_0", init_string)
  
      massweighter = np.repeat([mol.mass(a) for a in range(data["n_atom"])], 3)**(-0.5)
      B_pi = []
      H_pi = []
      irrep_lbls = mol.irrep_labels()
      max_disp = (findifrec["stencil_size"] - 1) // 2
      e_per_diag = 2 * max_disp
  
      # Unlike in the gradient case, we have no symmetry transformations to worry about.
      # We get to the task directly: assembling the force constants in each irrep block.
      for h in range(data["n_irrep"]):
          salc_indices = data["salc_indices_pi"][h]
          if not salc_indices: continue
  
          n_salcs = len(salc_indices)
          E = np.zeros((len(salc_indices), e_per_diag))
  
          # Step One: Diagonals
          # For asymmetric irreps, the energy at a + disp is the same as at a - disp
          # Just reuse the - disp energy for the + disp energy
  
          for i, salc_index in enumerate(salc_indices):
              for j in range(1, max_disp + 1):
                  E[i, max_disp - j] = displacements[f"{salc_index}: {-j}"]["energy"]
                  k = -j if h else j  # Because of the +- displacement trick
                  E[i, max_disp + j - 1] = displacements[f"{salc_index}: {k}"]["energy"]
          # Now determine all diagonal force constants for this irrep.
          if findifrec["stencil_size"] == 3:
              diag_fcs = E[:, 0] + E[:, 1]
              diag_fcs -= 2 * ref_energy
              diag_fcs /= (findifrec["step"]["size"]**2)
          elif findifrec["stencil_size"] == 5:
              diag_fcs = -E[:, 0] + 16 * E[:, 1] + 16 * E[:, 2] - E[:, 3]
              diag_fcs -= 30 * ref_energy
              diag_fcs /= (12 * findifrec["step"]["size"]**2)
          H_irr = np.diag(diag_fcs)
  
          # TODO: It's a bit ugly to use the salc indices to grab the off-diagonals but the indices
          # within the irrep to grab the diagonals. Is there a better way to do this?
  
          # Step Two: Off-diagonals
          # We need off-diagonal energies, diagonal energies, AND the reference energy
          # Grabbing off-diagonal energies is a pain, so once we know our SALCs...
          # ...define offdiag_en to do that for us.
          for i, salc in enumerate(salc_indices):
              for j, salc2 in enumerate(salc_indices[:i]):
                  offdiag_en = lambda index: displacements["{l}: {}, {k}: {}".format(k=salc, l=salc2, *data["disps"]["off"][index])]["energy"]
                  if findifrec["stencil_size"] == 3:
                      fc = (+offdiag_en(0) + offdiag_en(1) + 2 * ref_energy - E[i][0] - E[i][1] - E[j][0] - E[j][1]) / (
                          2 * findifrec["step"]["size"]**2)
                  elif findifrec["stencil_size"] == 5:
                      fc = (-offdiag_en(0) - offdiag_en(1) + 9 * offdiag_en(2) - offdiag_en(3) - offdiag_en(4) +
                            9 * offdiag_en(5) - offdiag_en(6) - offdiag_en(7) + E[i][0] - 7 * E[i][1] - 7 * E[i][2] +
                            E[i][3] + E[j][0] - 7 * E[j][1] - 7 * E[j][2] + E[j][3] + 12 * ref_energy) / (
                                12 * findifrec["step"]["size"]**2)
                  H_irr[i, j] = fc
                  H_irr[j, i] = fc
  
          B_pi.append(data["salc_list"].matrix_irrep(h))
          H_pi.append(_process_hessian_symmetry_block(H_irr, B_pi[-1], massweighter, irrep_lbls[h], data["print_lvl"]))
  
      # All blocks of the Hessian are now constructed!
      return _process_hessian(H_pi, B_pi, massweighter, data["print_lvl"])
function_tokens: def assemble_hessian_from_energies ( findifrec , freq_irrep_only ) : # This *must* be a Psi molecule at present - CdSalcList generation panics otherwise mol = core . Molecule . from_schema ( findifrec [ "molecule" ] , nonphysical = True , verbose = 0 ) displacements = findifrec [ "displacements" ] ref_energy = findifrec [ "reference" ] [ "energy" ] def init_string ( data ) : max_label_len = str ( max ( [ len ( label ) for label in displacements ] ) ) out_str = "" for label , disp_data in displacements . items ( ) : out_str += ( "    {:" + max_label_len + "s} : {:20.10f}\n" ) . format ( label , disp_data [ "energy" ] ) return ( "  Computing second-derivative from energies using projected, \n" "  symmetry-adapted, cartesian coordinates.\n\n" "  {:d} energies passed in, including the reference geometry.\n" "    Using {:d}-point formula.\n" "    Energy without displacement: {:15.10f}\n" "    Check energies below for precision!\n{}" . format ( len ( displacements ) + 1 , findifrec [ "stencil_size" ] , ref_energy , out_str ) ) data = _initialize_findif ( mol , freq_irrep_only , "2_0" , init_string ) massweighter = np . repeat ( [ mol . mass ( a ) for a in range ( data [ "n_atom" ] ) ] , 3 ) ** ( - 0.5 ) B_pi = [ ] H_pi = [ ] irrep_lbls = mol . irrep_labels ( ) max_disp = ( findifrec [ "stencil_size" ] - 1 ) // 2 e_per_diag = 2 * max_disp # Unlike in the gradient case, we have no symmetry transformations to worry about. # We get to the task directly: assembling the force constants in each irrep block. for h in range ( data [ "n_irrep" ] ) : salc_indices = data [ "salc_indices_pi" ] [ h ] if not salc_indices : continue n_salcs = len ( salc_indices ) E = np . zeros ( ( len ( salc_indices ) , e_per_diag ) ) # Step One: Diagonals # For asymmetric irreps, the energy at a + disp is the same as at a - disp # Just reuse the - disp energy for the + disp energy for i , salc_index in enumerate ( salc_indices ) : for j in range ( 1 , max_disp + 1 ) : E [ i , max_disp - j ] = displacements [ f"{salc_index}: {-j}" ] [ "energy" ] k = - j if h else j # Because of the +- displacement trick E [ i , max_disp + j - 1 ] = displacements [ f"{salc_index}: {k}" ] [ "energy" ] # Now determine all diagonal force constants for this irrep. if findifrec [ "stencil_size" ] == 3 : diag_fcs = E [ : , 0 ] + E [ : , 1 ] diag_fcs -= 2 * ref_energy diag_fcs /= ( findifrec [ "step" ] [ "size" ] ** 2 ) elif findifrec [ "stencil_size" ] == 5 : diag_fcs = - E [ : , 0 ] + 16 * E [ : , 1 ] + 16 * E [ : , 2 ] - E [ : , 3 ] diag_fcs -= 30 * ref_energy diag_fcs /= ( 12 * findifrec [ "step" ] [ "size" ] ** 2 ) H_irr = np . diag ( diag_fcs ) # TODO: It's a bit ugly to use the salc indices to grab the off-diagonals but the indices # within the irrep to grab the diagonals. Is there a better way to do this? # Step Two: Off-diagonals # We need off-diagonal energies, diagonal energies, AND the reference energy # Grabbing off-diagonal energies is a pain, so once we know our SALCs... # ...define offdiag_en to do that for us. for i , salc in enumerate ( salc_indices ) : for j , salc2 in enumerate ( salc_indices [ : i ] ) : offdiag_en = lambda index : displacements [ "{l}: {}, {k}: {}" . format ( k = salc , l = salc2 , * data [ "disps" ] [ "off" ] [ index ] ) ] [ "energy" ] if findifrec [ "stencil_size" ] == 3 : fc = ( + offdiag_en ( 0 ) + offdiag_en ( 1 ) + 2 * ref_energy - E [ i ] [ 0 ] - E [ i ] [ 1 ] - E [ j ] [ 0 ] - E [ j ] [ 1 ] ) / ( 2 * findifrec [ "step" ] [ "size" ] ** 2 ) elif findifrec [ "stencil_size" ] == 5 : fc = ( - offdiag_en ( 0 ) - offdiag_en ( 1 ) + 9 * offdiag_en ( 2 ) - offdiag_en ( 3 ) - offdiag_en ( 4 ) + 9 * offdiag_en ( 5 ) - offdiag_en ( 6 ) - offdiag_en ( 7 ) + E [ i ] [ 0 ] - 7 * E [ i ] [ 1 ] - 7 * E [ i ] [ 2 ] + E [ i ] [ 3 ] + E [ j ] [ 0 ] - 7 * E [ j ] [ 1 ] - 7 * E [ j ] [ 2 ] + E [ j ] [ 3 ] + 12 * ref_energy ) / ( 12 * findifrec [ "step" ] [ "size" ] ** 2 ) H_irr [ i , j ] = fc H_irr [ j , i ] = fc B_pi . append ( data [ "salc_list" ] . matrix_irrep ( h ) ) H_pi . append ( _process_hessian_symmetry_block ( H_irr , B_pi [ - 1 ] , massweighter , irrep_lbls [ h ] , data [ "print_lvl" ] ) ) # All blocks of the Hessian are now constructed! return _process_hessian ( H_pi , B_pi , massweighter , data [ "print_lvl" ] )