Geometry Optimisation

Geometry optimisation in ONETEP

This tutorial aims at showing how to run a simple geometry optimization with ONETEP. Geometry optimization is one of the primary tasks in quantum simulation. The essence of the calculation is for the constituting atoms to be moved to the positions where the total energy is minimal. In general, this can be tackled efficiently if the forces on the atoms can be computed.

Over the past twenty years, various schemes have been derived to solve this problem in the framework of ab initio calculations. These range from simple approaches based on molecular dynamics, such as the steepest descent and damped dynamics methods, to the more sophisticated conjugated gradient and Quasi-Newton methods.

The geometry optimization scheme implemented in ONETEP relies on the isolation of the atomic and electronic subsystems (i.e. the Born-Oppenheimer approximation). For a given configuration of the ionic positions, the electronic degrees of freedom are completely relaxed so that the electronic subsystem stays on the Born-Oppenheimer surface. All the possible configurations of the ionic positions therefore define a multi-dimensional potential energy surface for which we want to find the global minimum. The atomic forces are calculated by application of the Hellmann-Feynman theorem and the ionic positions are moved around by means of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method in order to find the minimum of the potential energy. At this point, one has to keep in mind that several local minima may be present in the configuration space and the algorithm can get trapped in one of those. Therefore, despite the sophistication of the minimization method, the location of a global minimum still requires the intuition of a good starting point.

The calculation flow of a geometry optimization in ONETEP is a three step process :

• Given an ionic configuration, the electronic degrees of freedom are relaxed (cfr. self-consistent optimization of the density kernel and NGWFs).
• The total energy and atomic forces are computed and compared with those of previous ionic configurations. The threshold chosen as stopping criterion for the geometry optimization is tested.
• The atomic position are updated by means of the BFGS algorithm.

The ethene molecule

As a first example of geometry optimisation, we will investigate the equilibrium configuration of the ethene molecule. You might consider to work in a new directory to run the following calculations.

As stated previously, the code relies on the pseudopotential approximation to describe the interaction between an ion and the electron cloud. You will need to copy the pseudo-potentials files corresponding to the hydrogen and carbon ions into the working directory.

Initially, we will rely on the default settings for the relaxation of the geometry. Therefore, the input file will be very similar to those created in the first tutorial. The only thing you will need is to activate the Geometry optimization scheme by setting task : GeometryOptimization. Try running a calculation with an energy cutoff of about 650 eV, NGWF radii of about 6.0 a0 and a cubic simulation cell of side-length 40 a0. For that purpose, use the keywords : cutoff_energy, lattice_cart, species, and species_pot.

A reasonable starting configuration for the ethene molecule is given here.

The calculation should take about 5min on 4 cpus. In the meantime you may want to repeat the procedure with varying parameters in order to converge the calculation with respect to the cutoff energy, the NGWF radii, as well as the size of the simulation cell. Besides, if you aim to compute a properties (e.g. the C-C bond length) with a given computational accuracy (e.g. 0.005 Ang), you should also check that the geom_max_iter and ngwf_threshold_orig parameters do not prevent to reach the desired accuracy.

The output of ONETEP consists principally of two files: ethene.out (the main output file) and ethene.geom. The main informations regarding the geometry optimization are gathered in ethene.geom. It contains one block of information for each iteration of the geometry optimization. Each block looks like :

 1 -1.36533889E+001 -1.36533889E+001 <-- E 4.00000000E+001 0.00000000E+000 0.00000000E+000 <-- h 0.00000000E+000 4.00000000E+001 0.00000000E+000 <-- h 0.00000000E+000 0.00000000E+000 4.00000000E+001 <-- h C 1 2.13562968E+001 1.98933496E+001 2.00755898E+001 <-- R H 1 2.23061517E+001 2.13601794E+001 2.10755433E+001 <-- R H 1 2.23049009E+001 1.84641257E+001 1.90930089E+001 <-- R C 1 1.86435063E+001 2.01062996E+001 1.99242388E+001 <-- R H 1 1.76950937E+001 2.15357977E+001 2.09069482E+001 <-- R H 1 1.76940505E+001 1.86382480E+001 1.89246710E+001 <-- R C 1 -1.50317812E-001 3.19801012E-002 -2.58367277E-004 <-- F H 1 1.68175906E-002 1.25101495E-002 1.04966796E-002 <-- F H 1 1.95828946E-002 -2.64393342E-002 -1.63619133E-002 <-- F C 1 1.50232461E-001 -3.20723476E-002 2.09285354E-004 <-- F H 1 -1.95654170E-002 2.63684588E-002 1.63201051E-002 <-- F H 1 -1.67497176E-002 -1.23470276E-002 -1.04057895E-002 <-- F

The first line is the iteration number.
The second line is the total energy.
The next three lines are the lattice vectors.
The next N lines (where N is the number of atoms) give the atomic coordinates.
The following N lines give the atomic forces.
All values are in Hartree atomic units.

You can use the perl script geom2xyz to generate a .xyz file containing the atomic coordinates at each iteration of the geometry optimization. Though the film of the relaxation provides you with crucial information such as the appearance of dissociation, symmetry breaking, etc..... It is a good practice to keep track of the energy and forces at each iteration in order to assess the relaxation process.

grep " E" ethene.geom | awk '{print \$1}' > ethene_energy.dat
gnuplot plot with lines 'ethene_energy.dat'

These commands should create a new file ethene_energy.dat and plot the evolution of the total energy. As expected, the total energy of the system decreases monotonically during the geometry optimization. Similaly the maximum rms force on the unconstrained atoms should decrease and reach the convergence threshold defined by geom_force_tol.

grep '|F|' ethene.out | awk '{print \$4}' > ethene_force.dat
gnuplot plot with lines 'ethene_force.dat'

Tuning the optimization procedure

You are now familiar with the geometry optimization scheme in ONETEP. You might examine in more details the input variables that allow to control the process. The keywords associated with the geometry optimization all start with the geom_ prefix. Take a few minutes to have a look at the variables : geom_max_iter, geom_convergence_win, geom_disp_tol, geom_energy_tol, and geom_force_tol.

Though their default values may appear to be convenient in most circumstances, these latter are the very basic input variables to master before launching a geometry optimization. It is important to note that the three tolerance criteria (geom_disp_tol, geom_energy_tol, and, geom_force_tol) are not exclusive. The three criteria have to be satisfied in order for the optimization to stop. You might have noticed that during the relaxation of ethene molecule, the default threshold imposed on the atomic forces (geom_force_tol : 0.02 Ha/Bohr) has been reached before the one associated with the convergence of the energy (geom_energy_tol : 10e- 06 Ha/Atom).

Like all the quasi-Newton schemes, the BFGS algorithm accumulates information about the Hessian matrix. As the the number of iteration increases, BFGS improves its knowledge of the potential energy surface around the minimum and the matrix used to build the quadratic model of the potential energy surface converges towards the true Hessian matrix at the local minimum. However, the Hessian matrix is poorly approximated during the first few relaxation steps. It is therefore important to properly initialize the BFGS scheme. This may be conveniently done by means of a unique parameter geom_frequency_est. For the best efficiency, its value should corresponds to the average of the optical phonon frequencies at the center of the Brillouin zone. In the case of the ethene molecule, the average of the experimentally reported vibration frequencies is 0.0081 Hartree. This value is very close to the default setting of geom_frequency_est and we do not expect any speed-up of the relaxation process by adjusting it.

In various circumstances, it may appears convenient to impose some constraints to the atomic positions during the geometry optimization. Note that in the case of molecular systems it is often a good idea to keep an atom or an axis fixed during the optimization process in order to avoid losing computational time due to the rotations and/or translations of the system. This can be done using the species_constraints keyword.

When running a geometry optimization, ONETEP produces a .continuation file. This latter contains all the information regarding the optimization process and can be very helpful to restart an optimization from a previous run. In such a case, the only thing you will need to do is to turn on the flag geom_continuation. In the same line of thought, you may want to avoid loosing time in writing the density kernels and NGWFs on the disk. An appropriate usage of the keywords write_converged_dkngwfs, read_denskern and read_tightbox_ngwfs may therefore help you to save some precious computational time.

More examples

At this point, you should be familiar with the keywords needed to run a proper geometry optimization. Therefore, we suggest you to leave the ethene molecule and to try to optimize a larger organic molecule. You will find an example input file for the sucrose molecule here

One notes that the write_converged_dkngwfs flag has been activated. In addition, the values of ngwf_cg_max_step and lnv_cg_max_step have been increased in order to allow unconstrained line search during the conjugate gradient optimization of the density kernel and NGWFs respectively. The calculation will take a while so is best launched on 8 or more cpus.

You should notice a rapid decrease of the maximum rms force on the ions during the first relaxation steps. However, the hydrogen atoms tend to wiggle quite a lot and it takes a some time for the positions to settle down according to relaxation criteria.

Here above, the geometry optimization scheme has been illustrated by means of two molecular systems. Obviously, the same scheme holds for periodic crystals. At last, we will investigate the adsorption of ammonia on a (10,8) carbon nanotube. You can find an example file containing the atomic positions of the nanotube as well as the adsorbate ammonia here. The carbon nanotube contains 488 carbon atoms in its unit-cell and its chiral periodicity is of 62.87 Bohr. Following the prescriptions stated above, you should be able to write an input file for the nanotube (an example file is given here).

Note that for large systems, the spatial extension of the density kernel has to be truncated in order to achieve the linear scaling. This can be done with the kernel_cutoff variable. Obviously, stringent truncation of the density kernel is expected to affect the accuracy of the calculation. Therefore, the cutoff length has to be carefully adjusted. Though the actual time needed for the relaxation procedure (on 16 cpus) will be of ~12 hours. You should already notice a significant decrease of the forces after the first few iterations.