Dear CP2K developers,<div> </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">CP2K paper</a> 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: ``` <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; overflow-wrap: break-word; border: 0px; border-radius: 0px; white-space-collapse: preserve; vertical-align: baseline;">[ 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: ``` <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; overflow-wrap: break-word; border: 0px; border-radius: 0px; white-space-collapse: preserve; vertical-align: baseline;">[ 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. </div> -- 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 <a href="mailto:cp2k+unsubscribe@googlegroups.com">cp2k+unsubscribe@googlegroups.com</a>. 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">https://groups.google.com/d/msgid/cp2k/285a20ce-76cd-454f-9441-ab6241e99c49n%40googlegroups.com</a>.