Problem with the calculation of the CUBE file and Wannier Centres using the ROKS method?

[Ilya Fedorov]

Dear colleagues,

 Could you please help me to understand the CP2K implementation of ROKS and the Wannier Centers.

I am facing the following two problems. I can't find an answer in the documentation. And, unfortunately, I don't really understand in the code.

1) The ROKS implementation in CP2K uses a spin-restricted calculation with multiple density matrices. As I understand it, this leads to some problems with printing CUBE files (for example, for each orbital it prints two cube files for spin_1 and spin_2). Which one is correct? Maybe the sum of them gives a correct density? 

(in log it print “Unclear how we define MOs in the restricted case ... skipping”)

 

2) I’m trying to use Wannier Centres with ROKS, but it doesn’t work in CP2K. As I understand this is also a multiple density matrix problem. I need Wannier Centres only (Many-body Position operator, no Jacobi rotation). Right now I use Wannier + ROKS in CPMD, and in CPMD the rotation is just switched off for the last two orbitals (diag(1,1)). 

Best regards

Ilya

[Ilya Fedorov]

I solved the problem and it was accepted into the code.

1. As I understand it, for ROKS the UKS base is used, and different spins of the same orbital have the same density. This allows to rely on the UKS codes, while you can use MO_CUBES only for spin-1. (regtest: QS/regtest-lsroks/O2_mo_cubes.inp)
   
   
2. I added such an implementation in the code. Only two method supports:
   • NONE: No rotation, just print the orbital centers.
   • JACOBI: Only doubly occupied orbitals rotate; for singly occupied orbitals, no rotation occurs (like the NONE method). (regtest: QS/regtest-lsroks/O2_loc_wan_jac.inp)

It is possible to use only NONE and JACOBI. Other methods can be added in the future, but ROKS calculates quite slowly, and speeding up the Wannier calculation does not make a big contribution here.

Commits: https://github.com/cp2k/cp2k/pull/3108

четверг, 24 августа 2023 г. в 15:19:45 UTC+3, Ilya Fedorov:

[Krack Matthias]

Hi Ilya

 

Thanks for contributing and fixing the print issue.

 

Indeed, ROKS produces only one set of orbitals, since there is also only one Hamiltonian, i.e. Kohn-Sham matrix, but different orbital occupations are applied for spin-up (alpha) and spin-down (beta) electrons which results in an alpha and beta electronic density matrix.

 

The LOW-SPIN ROKS implementation in CP2K based on OT plus ROTATION is experimental as indicated by a warning printed in the output. That implementation has not been maintained since a long time. Therefore, I am wondering, if the implementation works properly, e.g. for the new test case O2_mo_cubes.inp. I tried tighter thresholds and a larger cutoff and box size, but the SCF did not converge. What’s your experience and how do the results compare to ROKS in CPMD?

 

Best

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Ilya Fedorov <ilyafe...@gmail.com>
Date: Thursday, 16 November 2023 at 10:59
To: cp2k <cp...@googlegroups.com>
Subject: [CP2K:19518] Re: Problem with the calculation of the CUBE file and Wannier Centres using the ROKS method?

I solved the problem and it was accepted into the code.

1.       As I understand it, for ROKS the UKS base is used, and different spins of the same orbital have the same density. This allows to rely on the UKS codes, while you can use MO_CUBES only for spin-1. (regtest: QS/regtest-lsroks/O2_mo_cubes.inp)

2.       I added such an implementation in the code. Only two method supports:

o    NONE: No rotation, just print the orbital centers.

o    JACOBI: Only doubly occupied orbitals rotate; for singly occupied orbitals, no rotation occurs (like the NONE method). (regtest: QS/regtest-lsroks/O2_loc_wan_jac.inp)

It is possible to use only NONE and JACOBI. Other methods can be added in the future, but ROKS calculates quite slowly, and speeding up the Wannier calculation does not make a big contribution here.

 

Commits: https://github.com/cp2k/cp2k/pull/3108

четверг, 24 августа 2023 г. в 15:19:45 UTC+3, Ilya Fedorov:

Dear colleagues,

 Could you please help me to understand the CP2K implementation of ROKS and the Wannier Centers.

I am facing the following two problems. I can't find an answer in the documentation. And, unfortunately, I don't really understand in the code.

1) The ROKS implementation in CP2K uses a spin-restricted calculation with multiple density matrices. As I understand it, this leads to some problems with printing CUBE files (for example, for each orbital it prints two cube files for spin_1 and spin_2). Which one is correct? Maybe the sum of them gives a correct density? 

(in log it print “Unclear how we define MOs in the restricted case ... skipping”)

 

2) I’m trying to use Wannier Centres with ROKS, but it doesn’t work in CP2K. As I understand this is also a multiple density matrix problem. I need Wannier Centres only (Many-body Position operator, no Jacobi rotation). Right now I use Wannier + ROKS in CPMD, and in CPMD the rotation is just switched off for the last two orbitals (diag(1,1)). 

Best regards

Ilya

--
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...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/c479ab3a-9ed0-4750-9997-94a532c92a84n%40googlegroups.com.

[Ilya Fedorov]

Dear Matthias,

I am sorry for the delay in answering your question.

Indeed, sometimes I observed poor convergence of OT with ROKS, but sometimes the convergence is much better than diagonalization in СPMD. But at the moment I do not have a systematic understanding of the roots of this problem.

I made various comparisons between CPMD and CP2K in terms of calculations with ROKS: comparing forces, trajectories, orbitals and Wannier centers, in general I got a fairly close agreement.

In the near future I plan to collect these results into one report and post it here.

Best

Ilya

четверг, 16 ноября 2023 г. в 15:31:03 UTC+3, Krack Matthias:

[Ilya Fedorov]

Dear Matthias,

I apologise for such a long reply. The problem turned out to be not so simple.

At this moment it seems that ROKS does not work properly in CP2K.

An oxygen molecule turned out to be a poor example to use with ROKS. I just copied it from an existing regtest for ROKS, but CPMD and CP2K gave the S1 energies lower than the S0 energies.

A better example would be a formaldehyde molecule or a hydrogen molecule.The advantages of formaldehyde are that it has SOMO-1 and SOMO-2 separated in space, which, as will be discussed below, allows for better manifestation of the differences between CP2K and CPMD could be traced. I used a formaldehyde example from the original ROKS paper [Frank, I., Hutter, J., Marx, D., & Parrinello, M. (1998). JPC, 108, 4060]. 

Here are the main observations: 

1. The excitation energy E(S1) - E(S0) in both CPMD and CP2K packages is about 3.5 eV, which is in agreement with the paper of Frank et al.
2. The Wannier centres of all orbitals are almost the same between CPMD and CP2K (see Picture 1), but their order is different. This is important for SOMO-1 and SOMO-2, since we have the rotations turned off for these orbitals and their centres should have coincided.
3. The most important difference is in the visualisation of the orbitals. In the attached pdf-file (S0_S1_CP2K_CPMD.pdf) I visualise the orbitals for the ground state S0 and of the excited state S1. The isosurfaces are shown in red (+0.07) and in blue (-0.07). You can see that for the ground state the orbitals from CPMD and from CP2K coincide. For the excited state CPMD shows some change of order (with respect to S0) and the general picture remains close (the resemblance HOMO ~ SOMO-1 takes place). But for the S1 state in CP2K all orbitals look different, except last two (SOMO-1 and SOMO-2), and the last two orbitals seem to have switched places (the resemblance HOMO ~ SOMO-2 takes place) (the same pattern as with Wannier centres.).

I attach the codes for calculating the S0 and S1 states in CPMD and CP2K for formaldehyde. 

• The standard cpmd2cube.x script (https://github.com/CPMD-code/Addons/tree/main/cpmd2cube) should be used to convert the CPMD orbitals to .cube files (cmd: ./cpmd2cube.x WAVEFUNCTION.*)
• To calculate Wannier centres in CPMD, a separate code (cpmd_s1_wannier.inp) should be used which takes 1 MD step, this is because a separate case of calculating Wannier centres for ROKS is only implemented in MD. At the same time, in CPMD the calculation of Wannier centres itself affects the final output of orbitals

Best regards,

Ilya

P.S.
Although the energies of SOMO-1 and SOMO-2 judging by dE ~3.5 eV should be significantly different and could not simply get mixed up. In addition, it is worth noting that despite the complete difference in orbitals 1-5 in CP2K and CPMD, the final Wannier centres are the same, so we can talk about the influence of the difference of the basis or something else? 

Could it be that the OT method does not handle single occupied orbitals correctly or the ROKS method in CP2K does not handle double occupied orbitals correctly? Because in the case of CPMD we see that the shape of the double occupied orbitals has not changed, the order has changed and one new SOMO-2 orbital has been added. In CP2K, the new single occupied orbitals look correct, but the double occupied ones have coordinately changed in ROKS.

 

вторник, 21 ноября 2023 г. в 13:27:51 UTC+3, Ilya Fedorov:

[Ilya Fedorov]

Dear colleagues,

I found a bug and I have not corrected all issues with the implementation of orbitals printing with ROKS. I have now added make_mo_eig(), which returns orbitals from OT space. I have limited it to nspin=1 and everything now works consistently with CPMD. (https://github.com/cp2k/cp2k/pull/3274)

If you do not mind, may I replace the regtest with oxygen with formaldehyde, since it is a simpler example and has a larger dE in ROKS?

Best regards,
Ilya

понедельник, 12 февраля 2024 г. в 19:09:44 UTC+3, Ilya Fedorov: