[CP2K-user] [CP2K:19465] Computing MO coefficients using CP2K outputs

[mshakiba.k...@gmail.com]

Dear CP2K developers,

Hi, For educational purposes, I wanted to show how CP2K generates the MO 
coefficients. I've selected benzene with a SZV basis and the convergence 
limit is `1.0e-8`. Then, I printed the overlap (S), Kohn-Sham (K), and 
density matrices (P) in my input (using `&AO_MATRICES` section) and parsed 
them into Python. I used the same procedure as in CP2K paper 
<https://www.sciencedirect.com/science/article/abs/pii/S0010465505000615> as 
follows (I have attached the image as well):
```
# Cholesky decomposition of the Overlap matrix
U = np.linalg.cholesky(S).T
# U^1/2 was also tested
# U = scipy.linalg.fractional_matrix_power(S, 0.5)
# Printing of the U^T * U to make sure the decomposition is correct
print(np.dot(U.T, U))
# Computing U^-1 and U^T -1
U_inv = np.linalg.inv(U)
UT_inv = np.linalg.inv(U.T)
# Computing K'
K_prime = np.linalg.multi_dot( [UT_inv, K, U_inv] )
eigenvalues, eigenvectors = np.linalg.eig(K_prime)
sorted_indices = np.argsort(eigenvalues)
c = eigenvectors[:,sorted_indices]
# Back transformation
first_eigenvector = np.dot(U_inv, c[0])
```
Surprisingly, when I compare the outputs, I don't get the same coefficients 
as is printed out by CP2K in the MOLog files while I'm completely using 
CP2K data output. I think that the diagonalizer are similar and should not 
be different or at least not a lot but I don't know what is missing here.

For example, for the first MO, the numpy diagonalizer gives:
```
[ *0.24013208* 0.22018022 -0.31897544 0.31487772 0.03667346 -0.06520435 
-0.17002725 -0.08236077 0.07229922 -0.0671756 -0.01027799 -0.00122538 
-0.01032485 0.02372487 -0.04522648 0.07152516 -0.02957845 -0.04862388 
0.21459305 -0.2126356 -0.17272129 -0.87616482 0.55365848 -0.32681116 
-0.14926625 -0.16854633 -0.42447934 0.59950996 0.07111907 -0.2591646 ]
```
while CP2K output is:
```
[ *0.24013208* -0.04335288 -0.03266057 -0.02886998 0.25193464 -0.05879551 
-0.03198254 0.01560468 0.25274078 0.02363462 -0.03533591 -0.04778031 
0.27746047 -0.02251021 0.03359148 0.0590378 0.25630252 0.0619227 0.01419578 
-0.0243776 0.27224724 0.03675219 0.05132162 0.0233777 0.04862524 0.04652153 
0.04919569 0.05200084 0.04633218 0.05580193]```
You can see that the first element of the first eigenvector is the same 
while others are not. I also tried using higher number of digits while 
printing out CP2K output and also `scipy` diagonalizers but all give 
similar results.

It is worth noting that the eigenvectors indices in numpy are stored as the 
second index:
```
# Example
first_coeffs = eigenvectors[:,0]
```

For the convergence I compute the commutator relation and when computing 
the norm of the matrix I get a good convergence results in agreement with 
CP2K output:
e = KPS-SPK
```
e = np.linalg.multi_dot([K,P,S])-np.linalg.multi_dot([S,P,K])
print(np.linalg.norm(e))
>>> 1.2473125843097798e-08
```

I highly appreciate your time and any suggestion to get the correct results.

Thanks in advance.

-- 
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+unsubscribe at googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/782f3676-af56-4e46-b3ec-27f6f18c03d2n%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20231105/da80663d/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot 2023-11-05 203454.png
Type: image/png
Size: 130625 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20231105/da80663d/attachment-0001.png>