Structure optimization

In the tutorial on how to calculate atomization energies, we calculated the atomization energy for \(\rm{H}_2\) 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

Note

You must run the atomization script first.

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.