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

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

Hi Maximilian,

Thank you very much for your response. In fact, the problem was just solved 
:). I was making a trivial mistake in printing the eigenvectors. The only 
part that was needed to be fixed was by adding a transpose to the `c` 
matrix:

`c = eigenvectors[:,sorted_indices].T`

And it properly reproduces the coefficients and eigenvalues printed out by 
CP2K. I should have been more careful about that before posting the 
question. :)

Thanks again. 
On Monday, November 6, 2023 at 12:14:08 AM UTC-5 Maximilian Franz-Xaver 
Dorfner wrote:

> Hi,
> although I am no CP2K developer, I may can contribute to your answer. I 
> cannot comment your programm, since I have not looked deeply into it.
>
> To my impression benzene is a bad choice, when it comes to reproducing the 
> mo coefficients, as it hast lots of degeneracies, which reflect the 
> symmetry. At degeneracies the standard eigensolver routines converge to 
> some eigenvectors, but at say a twofold degeneracy the linear combination 
> of to eigenvector is (analytically) again an eigenvector. This implies 
> that, without an additional contstraint, it is very unlikely to reproduce 
> all the coefficients, as the the CP2K solver may converge to a rotated 
> solution compared to yours. 
>
> Best regards,
> Maximilian Dorfner 
>
>
>
>
>
> mshakiba.k... at gmail.com <mshakiba.k... at gmail.com> schrieb am Mo., 6. Nov. 
> 2023, 03:28:
>
>> 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.
>>
>> 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+uns... at googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/cp2k/285a20ce-76cd-454f-9441-ab6241e99c49n%40googlegroups.com 
>> <https://groups.google.com/d/msgid/cp2k/285a20ce-76cd-454f-9441-ab6241e99c49n%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
>>

-- 
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/53be5655-f944-4da9-b14e-a5b8dd7f6c02n%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20231105/924fe978/attachment-0001.htm>