<div dir="auto">Hi,<div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">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. </div><div dir="auto"><br></div><div dir="auto">Best regards,</div><div dir="auto">Maximilian Dorfner </div><div dir="auto"><br></div><div dir="auto"><br></div><div dir="auto"><br></div><br><br><div class="gmail_quote" dir="auto"><div dir="ltr" class="gmail_attr"><a href="mailto:mshakiba.k...@gmail.com">mshakiba.k...@gmail.com</a> <<a href="mailto:mshakiba.kerman.iran@gmail.com">mshakiba.kerman.iran@gmail.com</a>> schrieb am Mo., 6. Nov. 2023, 03:28:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear CP2K developers,<div><br></div><div>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 <a href="https://www.sciencedirect.com/science/article/abs/pii/S0010465505000615" target="_blank" rel="noreferrer">CP2K paper</a> as follows (I have attached the image as well):<br>```<br># Cholesky decomposition of the Overlap matrix<br>U = np.linalg.cholesky(S).T<br># U^1/2 was also tested<br># U = scipy.linalg.fractional_matrix_power(S, 0.5)<br># Printing of the U^T * U to make sure the decomposition is correct<br>print(np.dot(U.T, U))<br># Computing U^-1 and U^T -1<br>U_inv = np.linalg.inv(U)<br>UT_inv = np.linalg.inv(U.T)<br># Computing K'<br>K_prime = np.linalg.multi_dot( [UT_inv, K, U_inv] )<br>eigenvalues, eigenvectors = np.linalg.eig(K_prime)<br>sorted_indices = np.argsort(eigenvalues)<br>c = eigenvectors[:,sorted_indices]<br># Back transformation<br>first_eigenvector = np.dot(U_inv, c[0])<br>```<br>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.<br><br>For example, for the first MO, the numpy diagonalizer gives:<br>```<br><span style="box-sizing:border-box;overflow:auto;font-family:monospace;display:block;padding:1px 0px;margin:0px;line-height:inherit;color:rgb(0,0,0);word-break:break-all;border:0px;border-radius:0px;vertical-align:baseline">[ <b>0.24013208</b> 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 ] </span><br>```<br>while CP2K output is:<br>```<br><span style="box-sizing:border-box;overflow:auto;font-family:monospace;display:block;padding:1px 0px;margin:0px;line-height:inherit;color:rgb(0,0,0);word-break:break-all;border:0px;border-radius:0px;vertical-align:baseline">[ <b>0.24013208</b> -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]</span>```<br>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.<br><br>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:<br>e = KPS-SPK<br>```<br>e = np.linalg.multi_dot([K,P,S])-np.linalg.multi_dot([S,P,K])<br>print(np.linalg.norm(e))<br>>>> <span style="color:rgb(0,0,0);font-family:monospace">1.2473125843097798e-08</span><br>```<br><br>I highly appreciate your time and any suggestion to get the correct results.<br><br>Thanks in advance.<br><br></div> <p></p> -- <br> You received this message because you are subscribed to the Google Groups "cp2k" group.<br> To unsubscribe from this group and stop receiving emails from it, send an email to <a href="mailto:cp2k+unsubscribe@googlegroups.com" target="_blank" rel="noreferrer">cp2k+unsubscribe@googlegroups.com</a>.<br> To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/cp2k/285a20ce-76cd-454f-9441-ab6241e99c49n%40googlegroups.com?utm_medium=email&utm_source=footer" target="_blank" rel="noreferrer">https://groups.google.com/d/msgid/cp2k/285a20ce-76cd-454f-9441-ab6241e99c49n%40googlegroups.com</a>.<br><br> </blockquote></div></div> <p></p> -- <br /> You received this message because you are subscribed to the Google Groups "cp2k" group.<br /> To unsubscribe from this group and stop receiving emails from it, send an email to <a href="mailto:cp2k+unsubscribe@googlegroups.com">cp2k+unsubscribe@googlegroups.com</a>.<br /> To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/cp2k/CADd%2B%2BKVP7Y9ejJi2S3nNNsx9SRN4pXh4Zoo8UHkK7XWJq9zZ1w%40mail.gmail.com?utm_medium=email&utm_source=footer">https://groups.google.com/d/msgid/cp2k/CADd%2B%2BKVP7Y9ejJi2S3nNNsx9SRN4pXh4Zoo8UHkK7XWJq9zZ1w%40mail.gmail.com</a>.<br />