Results

Geometry optimization

It is well-know that the Kohn-Sham theorem is valid only for the ground state, which results in the charge density functional being very sensitive to the deviations from the ground state. This, in particular, can change the gap size or the Raman spectrum without even a hint to a mistake in the calculations. That's why it is important to make sure in our scf-run of Quantum Espresso that we have a good starting point. The best way to be sure is to do optimization by hands: check the system total energy  convergence with ecutwfc (which is known to be monotonous function) then find a good number of k-points (which is not monotonous in general) and, finally, to find the total energy minimum with respect to alat by fitting this dependence with a parabola in the vicinity of the potential equilibrium. However, as the number of optimized parameters increases, this procedure quickly becomes very tedious and time consuming, if not to say impossible for let say 5 optimization parameters.

On the other hand, this problem is a typical minimization task, though a bit non trivial as it deals with non-differentiable cost function that dependent on integer parameter that is a number of k-points. Luckily, there are methods that can reliably deal with such type of cost functions. One of the best is the Differential Evolution method for the global optimization. This method is designed to find a global minimum and as such is destined to find the best starting point for the Kohn-Sham density functional.

I am happy to share that I have realized Differential Evolution optimization as a python code and a single shell script creating and running python script with the code. The code is now available from GitHub repository: https://github.com/vasilsaroka/QEskillbox

Please, have also a look at Release notes for QEskillbox v0.1.0 - Geometry optimization https://github.com/vasilsaroka/QEskillbox/releases/tag/v0.1.0

You are welcome to use it. Any feedback will be much appreciated.

Band structure calculation

Having found a good approximation to the density functional, we usually want to calculate the energy band structure. Below the band structure of bulk Si is presented. This is a result of pw.x run calculation, data post processing with bands.x and visualization using gnuplot:

Effortless forceless geometry optimization

The standard methods of geometry optimization in classical and quantum molecular dynamics can be roughly divided into non-dumped and dumped dynamics methods. The first group encompasses such methods as the steepest descent, cojugate-gradient or bfgs, while the second includes different variations of the molecular dynamics in a viscose medium. Both approaches, however, are derivative-based methods where the direction of the descent is defined based on generalized forces calculated as the negative gradient of the potential (with some additional assumptions depending on the particular method). The non-dumped dynamic methods start from an approximate configuration of the system that is very close to the system potential energy minimum. Then the position of the atoms are promoted into the new position defined based on forces. In the dumped dynamic methods, the initial configuration is chosen as a starting point. The forces are calculated and the system is allowed to evolve according to the equations of motion in a viscous medium over a period of time. According to the equation of motion the system converts the excess of the potential energy into kinetic energy. Then the dumping parameter dissipates the kinetic energy of the system so that the system eventually descends to the potential energy minimum. In both cases, one can imagine a iterative system transformation as a local  force-driven evolution from the initial configuration to the final one. 

In contrast, Differential Evolution (DE) method represents a direct search, derivative-free technique. It starts from several randomly chosen configurations that evenly span the whole parameter search space. Each such configuration is described by a vector in the parameter space and the set of these vectors forms a population. The evolution of each vector in the population is defined by the differences between the parameter space vectors. Although it ostensibly can be related to the "approximate gradients" or "generalized forces", it is important to note that for a given vector this difference is taken between two other randomly chosen vectors within the population. This means the direction of descent is defined not in the point of the given vector but in the point that is somewhere far away in the parameter space. This later fact makes the dynamical picture given for the above derivative-based methods impossible. Despite this lack of "physical" meaning and clarity, the DE process shows a self-organizating feature. The whole iterative procedure looks like the convergence of all configurations in the population to a single configuration. Thus, forceless dynamics can be imagined as depicted in the animation below. At each iteration the best configuration can be chosen as the one with the minimum energy. 

In the figure above, the DE optimization of the silicene is presented. The population vectors in the parameter space are presented as several configurations of the real space lattices (light gray semi-transparent spheres). The the configuration with minimal energy is highlighted in brown. For this configuration, bonds between the atoms in the lattice are shown explicitly.  Since during the optimization the lattice constant is also optimized, the parameter space projection to the real space is slightly different for each configuration in the population as presented by gray and black boxes for the two atoms in the unit cell of the silicene. The parameter space projection to the real space for the configuration having minimal energy within the population is depicted with orange and red boxes. The total parameter space projection into a real space is the union of all the boxes.