docs: def constrained_opt_run(cls, vasp_cmd, lattice_direction, initial_strain,
                              atom_relax=True, max_steps=20, algo="bfgs",
                              **vasp_job_kwargs):
          """
          Returns a generator of jobs for a constrained optimization run. Typical
          use case is when you want to approximate a biaxial strain situation,
          e.g., you apply a defined strain to a and b directions of the lattice,
          but allows the c-direction to relax.
  
          Some guidelines on the use of this method:
          i.  It is recommended you do not use the Auto kpoint generation. The
              grid generated via Auto may fluctuate with changes in lattice
              param, resulting in numerical noise.
          ii. Make sure your EDIFF/EDIFFG is properly set in your INCAR. The
              optimization relies on these values to determine convergence.
  
          Args:
              vasp_cmd (str): Command to run vasp as a list of args. For example,
                  if you are using mpirun, it can be something like
                  ["mpirun", "pvasp.5.2.11"]
              lattice_direction (str): Which direction to relax. Valid values are
                  "a", "b" or "c".
              initial_strain (float): An initial strain to be applied to the
                  lattice_direction. This can usually be estimated as the
                  negative of the strain applied in the other two directions.
                  E.g., if you apply a tensile strain of 0.05 to the a and b
                  directions, you can use -0.05 as a reasonable first guess for
                  initial strain.
              atom_relax (bool): Whether to relax atomic positions.
              max_steps (int): The maximum number of runs. Defaults to 20 (
                  highly unlikely that this limit is ever reached).
              algo (str): Algorithm to use to find minimum. Default is "bfgs",
                  which is fast, but can be sensitive to numerical noise
                  in energy calculations. The alternative is "bisection",
                  which is more robust but can be a bit slow. The code does fall
                  back on the bisection when bfgs gives a non-sensical result,
                  e.g., negative lattice params.
              \*\*vasp_job_kwargs: Passthrough kwargs to VaspJob. See
                  :class:`custodian.vasp.jobs.VaspJob`.
  
          Returns:
              Generator of jobs. At the end of the run, an "EOS.txt" is written
              which provides a quick look at the E vs lattice parameter.
          """
          nsw = 99 if atom_relax else 0
  
          incar = Incar.from_file("INCAR")
  
          # Set the energy convergence criteria as the EDIFFG (if present) or
          # 10 x EDIFF (which itself defaults to 1e-4 if not present).
          if incar.get("EDIFFG") and incar.get("EDIFFG") > 0:
              etol = incar["EDIFFG"]
          else:
              etol = incar.get("EDIFF", 1e-4) * 10
  
          if lattice_direction == "a":
              lattice_index = 0
          elif lattice_direction == "b":
              lattice_index = 1
          else:
              lattice_index = 2
  
          energies = {}
  
          for i in range(max_steps):
              if i == 0:
                  settings = [
                          {"dict": "INCAR",
                           "action": {"_set": {"ISIF": 2, "NSW": nsw}}}]
                  structure = Poscar.from_file("POSCAR").structure
                  x = structure.lattice.abc[lattice_index]
                  backup = True
              else:
                  backup = False
                  v = Vasprun("vasprun.xml")
                  structure = v.final_structure
                  energy = v.final_energy
                  lattice = structure.lattice
  
                  x = lattice.abc[lattice_index]
  
                  energies[x] = energy
  
                  if i == 1:
                      x *= (1 + initial_strain)
                  else:
                      # Sort the lattice parameter by energies.
                      min_x = min(energies.keys(), key=lambda e: energies[e])
                      sorted_x = sorted(energies.keys())
                      ind = sorted_x.index(min_x)
                      if ind == 0:
                          other = ind + 1
                      elif ind == len(sorted_x) - 1:
                          other = ind - 1
                      else:
                          other = ind + 1 \
                              if energies[sorted_x[ind + 1]] \
                              < energies[sorted_x[ind - 1]] \
                              else ind - 1
                      if abs(energies[min_x]
                             - energies[sorted_x[other]]) < etol:
                          logger.info("Stopping optimization! Final %s = %f"
                                      % (lattice_direction, min_x))
                          break
  
                      if ind == 0 and len(sorted_x) > 2:
                          # Lowest energy lies outside of range of lowest value.
                          # we decrease the lattice parameter in the next
                          # iteration to find a minimum. This applies only when
                          # there are at least 3 values.
                          x = sorted_x[0] - abs(sorted_x[1] - sorted_x[0])
                          logger.info("Lowest energy lies below bounds. "
                                      "Setting %s = %f." % (lattice_direction, x))
                      elif ind == len(sorted_x) - 1 and len(sorted_x) > 2:
                          # Lowest energy lies outside of range of highest value.
                          # we increase the lattice parameter in the next
                          # iteration to find a minimum. This applies only when
                          # there are at least 3 values.
                          x = sorted_x[-1] + abs(sorted_x[-1] - sorted_x[-2])
                          logger.info("Lowest energy lies above bounds. "
                                      "Setting %s = %f." % (lattice_direction, x))
                      else:
                          if algo.lower() == "bfgs" and len(sorted_x) >= 4:
                              try:
                                  # If there are more than 4 data points, we will
                                  # do a quadratic fit to accelerate convergence.
                                  x1 = list(energies.keys())
                                  y1 = [energies[j] for j in x1]
                                  z1 = np.polyfit(x1, y1, 2)
                                  pp = np.poly1d(z1)
                                  from scipy.optimize import minimize
                                  result = minimize(
                                      pp, min_x,
                                      bounds=[(sorted_x[0], sorted_x[-1])])
                                  if (not result.success) or result.x[0] < 0:
                                      raise ValueError(
                                          "Negative lattice constant!")
                                  x = result.x[0]
                                  logger.info("BFGS minimized %s = %f."
                                              % (lattice_direction, x))
                              except ValueError as ex:
                                  # Fall back on bisection algo if the bfgs fails.
                                  logger.info(str(ex))
                                  x = (min_x + sorted_x[other]) / 2
                                  logger.info("Falling back on bisection %s = %f."
                                              % (lattice_direction, x))
                          else:
                              x = (min_x + sorted_x[other]) / 2
                              logger.info("Bisection %s = %f."
                                          % (lattice_direction, x))
  
                  lattice = lattice.matrix
                  lattice[lattice_index] = lattice[lattice_index] / \
                      np.linalg.norm(lattice[lattice_index]) * x
  
                  s = Structure(lattice, structure.species, structure.frac_coords)
                  fname = "POSCAR.%f" % x
                  s.to(filename=fname)
  
                  incar_update = {"ISTART": 1, "NSW": nsw, "ISIF": 2}
  
                  settings = [
                      {"dict": "INCAR",
                       "action": {"_set": incar_update}},
                      {"file": fname,
                       "action": {"_file_copy": {"dest": "POSCAR"}}}]
  
              logger.info("Generating job = %d with parameter %f!" % (i + 1, x))
              yield VaspJob(vasp_cmd, final=False, backup=backup,
                            suffix=".static.%f" % x,
                            settings_override=settings, **vasp_job_kwargs)
  
          with open("EOS.txt", "wt") as f:
              f.write("# %s energy\n" % lattice_direction)
              for k in sorted(energies.keys()):
                  f.write("%f %f\n" % (k, energies[k]))