Example: Epitaxial relaxation method

Absent dislocations, a film grown epitaxially on a substrate adopts the in-plane lattice parameters of that substrate. It is possible to model a sufficiently thick film as a bulk material on a virtual substrate. The trick is simply to lock the in-plane parameters while allowing the out-plane direction to relax.

VASP does not allow for this kind of relaxation directly. However, using Pylada when can simply design a method which will magically perform the relaxation. Below, a method based on the secant method is presented:

1. Find interval for which stress changes sign.
2. Bisect interval until required accuracy is obtained.

Admittedly, if one were smart, one could design a faster optimization method. Or one could use the optimization methods provided by scipy with a Vasp object, and be done with it, say this one.

First, we need a function capable of creating a strain matrix for the particular kind we have in mind, and apply it to the structure. The strain matrix can be obtained simply as the outer product of the epitaxial direction with itself.

def strain(x, direction, structure):
  """ Creates new structure with given input strain in c direction. """
  from numpy.linalg import inv
  from numpy import outer, dot
  result = structure.copy()
  # define strain matrix
  strain = outer(direction, direction) * x 
  # strain cell matrix and atoms
  result.cell = structure.cell + dot(strain, structure.cell)
  for atom in result: atom.pos += dot(strain, atom.pos)
  # return result.
  return result

The next step is to create a function which will compute the total energy of a strained structure while relaxing internal degrees of freedom. VASP is invoked in a separate sub-directory.

def function(vasp, structure, outdir, x, direction, **kwargs):
  from os.path import join
  directory = join(outdir, join("relax_ions", "{0:0<12.10}".format(x)))
  return vasp(strain(x, direction, structure), directory, relaxation='ionic', **kwargs)

The third component is a function to return the stress component in the epitaxial direction.

def component(extract, direction):
  """ Returns relevant stress component. """
  from numpy import dot
  return dot(dot(direction, extract.stress), direction)

Finally, we have all the pieces to put the bisecting algorithm together:

def epitaxial(vasp, structure, outdir=None, direction=[0,0,1], epiconv = 1e-4,
              initstep=0.05, **kwargs):
  # Compute initial structure.
  xstart = 0.0
  estart = function(vasp, structure, outdir, xstart, direction, **kwargs)
  
  # then checks stress for direction in which to release strain.
  stress_direction = 1.0 if component(estart, direction) > 0e0 else -1.0
  xend = initstep if stress_direction > 0e0 else -initstep
  eend = function(vasp, structure, outdir, xend, direction, **kwargs)

  # make sure xend is on other side of stress tensor sign.
  while stress_direction * component(eend, direction) > 0e0:
    xstart, estart = xend, eend
    xend += initstep if stress_direction > 0e0 else -initstep
    eend = function(vasp, structure, outdir, xend, direction, **kwargs)
  
  # now we have a bracket. We start bisecting it.
  while abs(estart.total_energy - eend.total_energy) > epiconv * float(len(structure)):
    # bisect interval and launch vasp.
    xmid = 0.5 * (xend + xstart)
    emid = function(vasp, structure, outdir, xmid, direction, **kwargs)
    # change interval depending on strain.
    if stress_direction * component(emid, direction) > 0: xstart, estart = xmid, emid
    else: xend, eend = xmid, emid

  # Finally, perform one last static calculation.
  efinal = eend if estart.total_energy > eend.total_energy else estart
  return vasp(efinal.structure, outdir=outdir, restart=efinal, relaxation='static', **kwargs)

Lines 4 through 16 correspond to (1) above. Starting from the initial input structure, we change the strain until an interval is found for which the stress component changes direction. The structure with minimum total energy will necessarily be found within this interval. Line 14 makes sure the interval is kept as small as possible.

We then move to the second part (lines 19-25) of the algorithm: bisecting the interval until the requisite convergence is achieved. Finally, one last static calculation is performed before returning. That’s probably an overkill, but it should be fairly cheap since we restart from the charge density and wavefunctions of a previous calculations.

The full listing of the code is given below. There are some extra bits to it, namely a docstring and some sanity checks on the function’s input parameter.

def strain(x, direction, structure):
  """ Creates new structure with given input strain in c direction. """
  from numpy.linalg import inv
  from numpy import outer, dot
  result = structure.copy()
  # define strain matrix
  strain = outer(direction, direction) * x 
  # strain cell matrix and atoms
  result.cell = structure.cell + dot(strain, structure.cell)
  for atom in result: atom.pos += dot(strain, atom.pos)
  # return result.
  return result

def function(vasp, structure, outdir, x, direction, **kwargs):
  from os.path import join
  directory = join(outdir, join("relax_ions", "{0:0<12.10}".format(x)))
  return vasp(strain(x, direction, structure), directory, relaxation='ionic', **kwargs)

def component(extract, direction):
  """ Returns relevant stress component. """
  from numpy import dot
  return dot(dot(direction, extract.stress), direction)

def epitaxial(vasp, structure, outdir=None, direction=[0,0,1], epiconv = 1e-4,
              initstep=0.05, **kwargs):
  """
      Performs epitaxial relaxation in given direction. 
  
      Performs a relaxation for an epitaxial structure on a virtual substrate.
      The external (cell) coordinates of the structure can only relax in the
      growth/epitaxial direction. Internal coordinates (ions), however, are
      allowed to relax in whatever direction. 
      
      Since VASP does not intrinsically allow for such a relaxation, it is
      performed by chaining different vasp calculations together. The
      minimization procedure itself is the secant method, enhanced by the
      knowledge of the stress tensor. The last calculation is static, for
      maximum accuracy.

      :param vasp: 
        :py:class:`Vasp <pylada.vasp.Vasp>` functional with wich to perform the
        relaxation.
      :param structure:
        :py:class:`Structure <pylada.crystal.Structure>` for which to perform the
        relaxation.
      :param str outdir: 
        Directory where to perform calculations. If None, defaults to current
        working directory. The intermediate calculations are stored in the
        relax_ions subdirectory.
      :param direction:
        Epitaxial direction. Defaults to [0, 0, 1].
      :param float epiconv: 
        Convergence criteria of the total energy.
  """
  from os import getcwd
  from copy import deepcopy
  from numpy.linalg import norm
  from numpy import array

  # takes care of input parameters.
  if outdir is None: outdir = getcwd
  vasp = deepcopy(vasp) # keep function stateless!
  direction = array(direction) / norm(direction)
  kwargs.pop('relaxation', None) # we will be saying what to optimize.
  if 'ediff' in kwargs: vasp.ediff = kwargs.pop('ediff')
  if vasp.ediff < epiconv: vasp.ediff = epiconv * 1e-2

  # Compute initial structure.
  xstart = 0.0
  estart = function(vasp, structure, outdir, xstart, direction, **kwargs)
  
  # then checks stress for direction in which to release strain.
  stress_direction = 1.0 if component(estart, direction) > 0e0 else -1.0
  xend = initstep if stress_direction > 0e0 else -initstep
  eend = function(vasp, structure, outdir, xend, direction, **kwargs)

  # make sure xend is on other side of stress tensor sign.
  while stress_direction * component(eend, direction) > 0e0:
    xstart, estart = xend, eend
    xend += initstep if stress_direction > 0e0 else -initstep
    eend = function(vasp, structure, outdir, xend, direction, **kwargs)
  
  # now we have a bracket. We start bisecting it.
  while abs(estart.total_energy - eend.total_energy) > epiconv * float(len(structure)):
    # bisect interval and launch vasp.
    xmid = 0.5 * (xend + xstart)
    emid = function(vasp, structure, outdir, xmid, direction, **kwargs)
    # change interval depending on strain.
    if stress_direction * component(emid, direction) > 0: xstart, estart = xmid, emid
    else: xend, eend = xmid, emid

  # Finally, perform one last static calculation.
  efinal = eend if estart.total_energy > eend.total_energy else estart
  return vasp(efinal.structure, outdir=outdir, restart=efinal, relaxation='static', **kwargs)