Global optimization of correlation matrices for DFT+U (Abinit)¶
The variable dmatpawu¶
The objective of this global search is finding the optimal values for the density matrices in DFT+U. ABINIT allows to locally optimize the density matrices from a given initial value. The initial density matrices used in LDA+U are kept fixed during the first usedmatpu SCF iterations. For SCF iterations beyond that, the density matrices change. The challenge here is that from a given initial set of density matrices the system gets easily trapped into a local minimun, the usual procedure then is to start from several initial options hoping to reach the global minimum at some point.
A global minimizer takes the responsability of efficiently explore the configuration space of the problem. PyChemia implements several global searchers as we saw on previous tutorials. Those global searchers joined by an efficient evaluation infraestructure allows many evaluations being perform without human assistance and an effective chance of reaching a global minimum if the search is long enough.
For this tutorial consider the following problem:
‘pychemia/test/data/abinit_dmatpaw/abinit.in’.
First, we can read the abinit input and access the contents of the variable ‘dmatpawu’:
import pychemia
import numpy as np
pychemia_path = pychemia.__path__[0]
abiinput = pychemia.code.abinit.InputVariables(pychemia_path + '/test/data/abinit_dmatpawu/abinit.in')
dmatpawu = np.array(abiinput['dmatpawu']).reshape(-1,5,5)
The variable ‘dmatpawu’ stores the contents of 4 5x5 matrices, the correlation matrices for the corresponding 4 Co atoms in the crystal. The matrices can be converted into a numpy array with shape (4, 5, 5) and they look like this:
array([[[ 0.06256, 0. , 0.01218, 0. , 0. ],
[ 0. , 0.481 , 0. , 0.45877, 0. ],
[ 0.01218, 0. , 0.07148, 0. , 0. ],
[ 0. , 0.45877, 0. , 0.481 , 0. ],
[ 0. , 0. , 0. , 0. , 0.94038]],
[[ 0.97852, 0. , 0.00305, 0. , 0. ],
[ 0. , 0.95493, 0. , -0.00449, 0. ],
[ 0.00305, 0. , 0.98115, 0. , 0. ],
[ 0. , -0.00449, 0. , 0.95493, 0. ],
[ 0. , 0. , 0. , 0. , 0.95168]],
[[ 0.97852, 0. , -0.00305, 0. , 0. ],
[ 0. , 0.95493, 0. , 0.00449, 0. ],
[-0.00305, 0. , 0.98115, 0. , 0. ],
[ 0. , 0.00449, 0. , 0.95493, 0. ],
[ 0. , 0. , 0. , 0. , 0.95168]],
[[ 0.06256, 0. , -0.01218, 0. , 0. ],
[ 0. , 0.481 , 0. , -0.45877, 0. ],
[-0.01218, 0. , 0.07148, 0. , 0. ],
[ 0. , -0.45877, 0. , 0.481 , 0. ],
[ 0. , 0. , 0. , 0. , 0.94038]]])
The objective is to find the set of correlation matrices that minimize the energy. Those are density matrices so even if we have 100 numbers, any set of numbers is a valid set of correlation matrices. We will now convert this set of matrices into a reduced set of variables that can be treated independently.
A correlation matrix can be express as the following product:
R*O*R^{-1}
Where R is a rotation matrix, O is a diagonal matrix with a trace that is the total number of electrons correlated. We need to find a set of independent variables to recreate any correlation matrix. We know that not any arbitrary set of 25 numbers is a good rotation matrix. However, a 5x5 rotation matrix can be effectively decomposed into 10 independent numbers, the so called “Generalized Euler angles”, this set of angles reduces the 25 values from a 5x5 rotation matrix into 10 independent variables. We should also be aware that the occupations on the diagonal of the matrix ‘O’ are not exactly integers, we will account for the small differences into a separate set of values. With those premises a 5x5 correlation matrix is converted into a set with 10 euler angles, 5 occupations and 5 deltas. This is done by the routines ‘dmatpawu2params’ and ‘params2dmatpawu’ that allow us to go back and forward from the set of correlation matrices into the set of ‘euler_angles’, intger ‘occupations’ and the small ‘deltas’:
params = pychemia.population.orbitaldftu.dmatpawu2params(dmatpawu, 5)
The variable ‘params’ is a dictionary with values for ‘occupations’, ‘deltas’ and ‘euler_angles’:
{'deltas': array([[ 0.02223 , 0.054049, 0.079991, 0.06023 , 0.05962 ],
[ 0.04956 , 0.04832 , 0.04058 , 0.023486, 0.016844],
[ 0.04956 , 0.04832 , 0.04058 , 0.023486, 0.016844],
[ 0.02223 , 0.054049, 0.079991, 0.06023 , 0.05962 ]]),
'euler_angles': array([[ 0. , 0. , 0. , 0. , 0.785398, 0. ,
0. , -0.609893, 0. , -1.570796],
[-0.581865, 0. , -1.570796, 0. , 0.785398, 0. ,
1.570796, 1.570796, 1.570796, 3.141593],
[-0.581865, 0. , 1.570796, 0. , 0.785398, 0. ,
1.570796, 1.570796, 1.570796, -0. ],
[ 0. , 0. , 0. , 0. , 0.785398, 0. ,
0. , -0.609893, 0. , 1.570796]]),
'ndim': 5,
'num_matrices': 4,
'occupations': array([[0, 0, 0, 1, 1],
[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[0, 0, 0, 1, 1]])}
This is in fact the set of independent variables that we can use to optimize the correlation using a global searcher. We can also go back and recover the correlation matrices using the inverse procedure, ie, from the dictionary params recover the correlation matrices:
dmatpawu_new = pychemia.population.orbitaldftu.params2dmatpawu(params)
The dmatpawu is recovered from the values stored in ‘params’:
array([[[ 0.06256, -0. , -0.01218, -0. , 0. ],
[-0. , 0.481 , -0. , 0.45877, 0. ],
[-0.01218, -0. , 0.07148, -0. , 0. ],
[-0. , 0.45877, -0. , 0.481 , 0. ],
[ 0. , 0. , 0. , 0. , 0.94038]],
[[ 0.97852, -0. , 0.00305, -0. , -0. ],
[-0. , 0.95493, -0. , -0.00449, -0. ],
[ 0.00305, -0. , 0.98115, -0. , -0. ],
[-0. , -0.00449, -0. , 0.95493, 0. ],
[-0. , -0. , -0. , 0. , 0.95168]],
[[ 0.97852, -0. , -0.00305, -0. , -0. ],
[-0. , 0.95493, 0. , 0.00449, 0. ],
[-0.00305, 0. , 0.98115, -0. , 0. ],
[-0. , 0.00449, -0. , 0.95493, 0. ],
[-0. , 0. , 0. , 0. , 0.95168]],
[[ 0.06256, 0. , -0.01218, -0. , 0. ],
[ 0. , 0.481 , 0. , -0.45877, 0. ],
[-0.01218, 0. , 0.07148, -0. , 0. ],
[-0. , -0.45877, -0. , 0.481 , 0. ],
[ 0. , 0. , 0. , 0. , 0.94038]]])
Each correlation matrix contains 25 values, using the procedure above, we reduce this number to 20: 10 euler angles, 5 integer occupations and 5 deltas. The values of deltas can be ignored for the purpose of the global searcher and the occupations are contrained by the condition that their sum must be the equal to the number of electrons in the correlated orbital. We have now the ingredients to move into the next step, create a population of correlation matrices.
The population¶
The most simple way of creating the population requires just the name of the mongo database to be created and one abinit input file. The relevant information to setup the search will be infered from the contents of the abinit input file:
popu=pychemia.population.orbitaldftu.OrbitalDFTU('test', abinit_input=pychemia_path +
'/test/data/abinit_dmatpawu/abinit.in')
Orbital population:
Species [znucl]: [19, 27, 9]
Orbitals corrected:
19 : False
27 : True (l=2)
9 : False
Number of atoms where DFT+U is applied: 4
Correlation of 'd' orbitals
Variables controling the total number of matrices
nsppol : 1
nspinor: 1
nspden : 2
Total number of matrices expected on dmatpawu: 4
Number of electrons for each correlation matrix: [2 5 5 2]
Number of independent matrices: 4
Create random correlation matrices can be done with:
popu.add_random()
The return is the Indentifier of the new entry on the database. Also a set of new random correlation matrices can be created with:
popu.random_population(16)
We have the basic ingredients for creating the first population for the global searcher. How the correlation matrices are evaluated is out of scope of the population and depends on the particularities of the machines where Abinit is used to evaluate them. We will move our focus to the methods needed to produced new correlation matrices based on the results of a given set of correlation matrices.