Computing the ground state of SiH4 using direct minimization

Instead of using the scf function to compute the ground state energy and electron density of a molecule or solid by via a self-consistent field iteration, we can also use the trdcm function to compute the ground state via a constrained minimization algorithm. This function is particular useful for systems with small gaps at zero temperature (See the subsection on regularization under the Algorithm section). The use of a trust region in the constrained minimization procedure stablizes the optimization algorithm so that no charge sloshing occurs. However, for systems with very small gaps, the convergence can be slow or the optimization procedure can stagnate due to the ill-posed nature of the problem.

Just like the sih4_scf example describe in the “Getting started” section, we set up the SiH4 molecule and its supercell by the following lines of codes

%
% construct the SiH4 (Silane) molecule
%
kssolvpptype('default');
%
% 1. construct atoms
%
a1 = Atom('Si');
a2 = Atom('H');
atomlist = [a1 a2 a2 a2 a2];
%
% 2. set up supercell
%
C = 10*eye(3);
%
% 3. define the coordinates the atoms
%
coefs = [
 0.0     0.0      0.0
 0.161   0.161    0.161
-0.161  -0.161    0.161
 0.161  -0.161   -0.161
-0.161   0.161   -0.161
];
xyzlist = coefs*C';
%
% 4. Configure the molecule (crystal)
%
mol = Molecule('supercell',C, 'atomlist',atomlist, 'xyzlist' ,xyzlist, ...
    'ecutwfc',12.5, 'name','SiH4' );

We then call the trdcm function with the default algorithmic settings by

[mol,H,X,info] = trdcm(mol);

The trdcm run converges in 10 iterations. At convergence, it prints out the following energy and error statistics

DCM iter = 10
Rel Vtot Err    =            1.708e-06
Total Energy    = -6.2288552012953e+00
Etot            = -6.2288552012953e+00
Ekin            = -5.3005625021617e+00
Eewald          = -1.5452143529686e+00
Ealphat         =  0.0000000000000e+00
Ecor            = -2.5585860020182e+00
Ecoul           =  3.2037896636954e+00
Exc             = -2.5868680098604e+00
--------------------------------------
Total time used =            1.298e+01
||HX-XD||_F     =            5.120e-06