Structure optimization

In this tutorial we consider structure optimization of the H₂ molecule. For that, we will calculate the atomization energy of the molecule both for the experimentally determined geometry and for the structure relaxed using the GPAW calculator.

In a previous tutorial we have calculated atomization energies. As a short recap, the following script will calculate the atomization energy of a hydrogen molecule with the geometry from experiment:

from ase import Atoms
from gpaw import GPAW

a = 10.  # Size of unit cell (Angstrom)
c = a / 2

# Hydrogen atom:
atom = Atoms('H',
             positions=[(c, c, c)],
             magmoms=[0],
             cell=(a, a + 0.0001, a + 0.0002))  # break cell symmetry

# gpaw calculator:
calc = GPAW(mode='pw',
            xc='PBE',
            hund=True,
            txt='H.out')
atom.calc = calc

e1 = atom.get_potential_energy()
calc.write('H.gpw')

# Hydrogen molecule:
d = 0.74  # experimental bond length
molecule = Atoms('H2',
                 positions=([c - d / 2, c, c],
                            [c + d / 2, c, c]),
                 cell=(a, a, a))

calc = calc.new(hund=False,  # no hund rule for molecules
                txt='H2.out')

molecule.calc = calc
e2 = molecule.get_potential_energy()
calc.write('H2.gpw')

Above, we calculated the atomization energy for H₂ using the experimental bond length of 0.74 Å. In this tutorial, we ask an ASE optimizer to iteratively find the structural energy minimum, where all atomic forces are below 0.05 eV/Å. The following script will do the job:

# web-page: optimization.txt

from gpaw import restart
from ase.parallel import paropen as open
from ase.optimize import QuasiNewton


molecule, calc = restart('H2.gpw', txt='H2-relaxed.txt')

e2 = molecule.get_potential_energy()
d0 = molecule.get_distance(0, 1)

with open('optimization.txt', 'w') as fd:
    print('experimental bond length:', file=fd)
    print(f'hydrogen molecule energy: {e2:5.2f} eV', file=fd)
    print(f'bondlength              : {d0:5.2f} Ang', file=fd)

    # Find the theoretical bond length:
    relax = QuasiNewton(molecule, logfile='qn.log')
    relax.run(fmax=0.05)

    e2 = molecule.get_potential_energy()
    d0 = molecule.get_distance(0, 1)

    print(file=fd)
    print('PBE energy minimum:', file=fd)
    print(f'hydrogen molecule energy: {e2:5.2f} eV', file=fd)
    print(f'bondlength              : {d0:5.2f} Ang', file=fd)

The result is:

experimental bond length:
hydrogen molecule energy: -6.65 eV
bondlength              :  0.74 Ang

PBE energy minimum:
hydrogen molecule energy: -6.65 eV
bondlength              :  0.75 Ang

To save time you could have told the minimizer to keep one atom fixed, and only relaxing the other. This is achieved through the use of constraints:

molecule.set_constraint(FixAtoms(mask=[0, 1]))

The keyword mask contains list of booleans for each atom indicating whether the atom’s position should be fixed or not. See the ase.constraints module on the ASE page for more information and examples for setting constraints.