Getting started¶
This section shows how one can use KSSOLV2.0 to set up and solve a simple problem.
To get started, let us start a MATLAB interactive session.
Move into the directory in which KSSOLV2.0 is installed. The directory should contain the file
KSSOLV_Startup.m
Type ‘KSSOLV_Startup’ (without the single quotes) at the MATLAB command window:
>> KSSOLV_Startup
This m-file makes all KSSOLV2.0 functions in different directories accessible from the command line.
You can then go into the example directory, and run the sih4_scf.m example.
>> sih4_scf
This simple scripts performs a self-consistent field (SCF) iteration to find the ground state energy and electron density of a silane molecule.
The output of the run looks like the following:
The pseudopotential for H is loaded from C:\cygwin64\home\CYang\kssolv_sandbox/ppdata/default/H_ONCV_PBE-1.0.upf
The pseudopotential for Si is loaded from C:\cygwin64\home\CYang\kssolv_sandbox/ppdata/default/Si_ONCV_PBE-1.0.upf
Regular SCF for Pure DFT
Beging SCF calculation for SiH4...
SCF iter 1:
Rel Vtot Err = 9.581e-02
Total Energy = -6.2142942363861e+00
SCF iter 2:
Rel Vtot Err = 1.850e-02
Total Energy = -6.2284207490181e+00
...
SCF iter 9:
Rel Vtot Err = 7.816e-09
Total Energy = -6.2288552013125e+00
Convergence is reached!
Elapsed time is 1.627365 seconds.
Etot = -6.2288552013125e+00
Eone-electron = -5.3005632795888e+00
Ehartree = 3.2037906874381e+00
Exc = -2.5868682561932e+00
Eewald = -1.5452143529686e+00
Ealphat = 0.0000000000000e+00
--------------------------------------
Total time used = 5.906e+00
||HX-XD||_F = 7.884e-09
Let’s now take a closer look at the script sih4_scf.m to see how we can easily set up a molecule and perform an SCF calculation to find its ground state.
In the first few lines of the code, we define both the silicon and hydrogen atoms and put one silicon and four hydrogen atoms in an array of atoms named atomlist.
%
% construct the SiH4 (Silane) molecule
%
% 1. construct atoms
%
a1 = Atom('Si');
a2 = Atom('H');
atomlist = [a1 a2 a2 a2 a2];
We then define a few things that are required to construct the silane molecule.
%
% 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';
The array C defines the 10x10x10 (Bohr^3) supercell that contains the silane molecule. The coefs array gives the position of each atom relative to the size of the supercell. These relative positions are converted to absolute coordinates (in Bohr) through the matrix-matrix multiplication coefs*C’, which is the same as (C*coefs'). Note that the order of the coordinate list in coefs must be consistent with the order of the atoms in the atomlist array. Also, the coordinates can be negative because the supercell is assumed to be periodically extended in a planewave calculation.
The arrays C, atomlist and xyzlist are used to define a
Molecule object named mol in the following segment of the script.
%
% 4. Configure the molecule (crystal)
%
mol = Molecule('supercell',C, 'atomlist',atomlist, 'xyzlist' ,xyzlist, ...
'ecutwfc',12.5, 'name','SiH4' );
Also included in the attributes of the Molecule object are the planewave cutoff energy ecutwfc, which is set to 12.5 Hatree and an optional name: SiH4.
The planewave cutoff energy is a discretization parameter that determines the accuracy of the approximation and the cost of the computation. The larger the cutoff energy, the more accurate the approximation of the DFT energy to the true ground state energy of the molecule, and the higher the cost of the computation.
Before we perform an SCF calculation, we set a few optional parameters by using the setksopt function.
options = setksopt();
options.maxscfiter = 20;
%options.verbose = 1;
The setksopt() call returns a structure options with a set of default parameters used by the SCF calculation. We can reset them by either calling the function again with appropriate (key,value) pairs, or simply reset a particular field of the structure. e.g., options.maxiscfiter=20 resets the the maximum number of SCF iterations allow to 20 from the default value of 10. We can also reset options.verbose to 1 to see more intermediate output from the SCF run.
Finally, we can call the scf function with the constructed Molecule object named mol and option strcuture options by
[mol,H,X,info] = scf(mol,options);
The function returns an updated Molecule object mol, a Kohn-Sham Hamiltonian object H corresponding to the electron density and total potential produced in the last SCF iterations, a Kohn-Sham wavefunction object X which contains Kohn-Sham orbitals (or coefficients of the planewave expansion of the orbitals), and some additional information contained in the info structure.
More information about H and X can be found in the Class section of the documentation.