Structure relaxation recipes

In the previous sections we reviewed the basics of structure relaxations. Here, we give some examples how to run and refine structure optimizations with GPAW.

Fixed-cell relaxation

The simplest scenario is to optimize the atom positions without changing the unit cell shape or size. In this case, we need to calculate the forces acting on the atoms and pass these to an optimization routine (here: BFGS) which iteratively updates the atom positions such that the magnitude of the forces decreases jointly with the total potential energy. First import the necessary modules:

from ase.build import bulk
from ase.optimize import BFGS
from gpaw import GPAW

Now consider the function relax:

• expects an ASE Atoms object, and

• calculation parameters dictionary param,

• optionally the maximum force fmax can be set,

• extracts calculation parameters from param (selected by the keyword fast_forces),

• sets up a GPAW calculator,

• adds the initial magnetic moments,

• optionally uses van der Waals forces with the DFT-D3 method, and

• sets up the optimization.

def relax(atoms, calculator_params,
          fmax=0.01, d3=False, fixcell=True,
          logname='opt.log',
          trajname='opt.traj'):

    # set DFT calculator
    calc_dft = GPAW(**calculator_params)

    # magnetize atoms
    atoms.set_initial_magnetic_moments(len(atoms) * [1])
    # non-magnetic calculation:
    # atoms.set_initial_magnetic_moments(len(atoms) * [0])

    # optionally include van der Waals DFT-D3
    if d3:
        from ase.calculators.dftd3 import DFTD3
        calc = DFTD3(dft=calc_dft)
    else:
        calc = calc_dft

    # set calculator
    atoms.calc = calc

    # set configuration to be optimized
    if fixcell:
        # only optimize positions of the atoms
        opt_conf = atoms
    else:
        # setup full relaxation
        # set unit cell filter
        opt_conf = FrechetCellFilter(atoms)

    # setup optimizer
    # specify logfile and trajectory file names
    opt = BFGS(opt_conf, logfile=logname, trajectory=trajname)
    # run the optimization until forces are smaller than fmax
    opt.run(fmax=fmax)

    return atoms

The relaxation history is written to logname = 'opt.log', its trajectory is saved to trajname = 'opt.traj' and the final configuration is returned from the function.

The parameter file params_forces.py contains the calculation parameters in standardized format. Typically we would start with a fast and less accurate relaxation, which is consecutively refined with tighter convergence parameters once we approach the atomic configuration which minimizes the potential energy (see accuracy of the self-consistency cycle and converging forces for relevant options).

# fast forces
calculator_params = {
    "xc": "PBE",
    "basis": "dzp",
    "mode": {"name": "lcao"},
    "kpts": {"size": [1, 1, 1],
             "gamma": True},
    "convergence": {"density": 1e-3,
                    "forces": 1e-2},
    "occupations": {"name": "fermi-dirac",
                    "width": 0.05},
    "mixer": {"method": "fullspin",
              "backend": "pulay"},
    "txt": "rlx.txt",
}

For more accurate force convergence the Brillouin zone sampling should be refined by setting a larger number of kpoints (e.g. five or more points in each direction), and the convergences criteria should be tightened to "convergence": {"density": 1e-6, "forces": 1e-4}.

Full relaxation

The cell shape and size can be relaxed together with the atom positions using cell filters. For that, we need to import the corresponding class using from ase.filters import FrechetCellFilter. The cell filter takes the Atoms object as input and is directly handed to the optimization routine. For full relaxations, you can then simply replace the optimization configuration line opt_conf = atoms with the opt_conf = FrechetCellFilter(atoms). In the example script we use an if-statement for that.

Note, that calculation of stresses requires the plane wave mode in GPAW and that the accuracy of the calculations should be refined (here: taking the more accurate parameter set params_stresses.py).

# accurate stress
calculator_params = {
    "xc": "PBE",
    "basis": "dzp",
    "mode": {"name": "pw",
             "ecut": 600},
    "kpts": {"size": [5, 5, 5],
             "gamma": True},
    "convergence": {"density": 1e-6,
                    "forces": 1e-4},
    "occupations": {"name": "fermi-dirac",
                    "width": 0.05},
    "mixer": {"method": "fullspin",
              "backend": "pulay"},
    "txt": "rlx.txt",
}

Here, you can download the full relaxation script relax.py.