DFT+U unphysical warning

724 views
Skip to first unread message

Nicholas Winner

unread,
Apr 16, 2021, 4:56:47 PM4/16/21
to cp2k

Hello all, I am running some DFT+U calculations on Mn-O systems. While I have found a cp2k effective U value in the literature of ~1.3 eV for Mn. I notice that when tuning the value myself, I begin seeing the following warning once the U value reaches 0.5eV.

*** WARNING in dft_plus_u.F:2006 :: DFT+U energy contibution is negative ***

 *** possibly due to unphysical Mulliken charges!    

Now this is only a warning, not an indication that the calculation is *necessarily* wrong, but it is troubling at least. Especially when my U value is nowhere near the size of lit value. I am using PLUS_U_METHOD MULLIKEN_CHARGES in order to have a marginally more robust solution. Does anyone have experience with how seriously to take this warning? I don't have a frame of reference to know if I should ignore it.


-Nick

Krack Matthias (PSI)

unread,
Apr 17, 2021, 2:00:09 AM4/17/21
to cp...@googlegroups.com

Hello Nick

 

I do not recommend the use of the PLUS_U_METHOD MULLIKEN_CHARGES. Just use the PLUS_U_METHOD MULLIKEN. The Mulliken population analysis can give unphysical orbital occupations, i.e. values greater than one (UKS case) or two (RKS case). Often the maximum occupation is only slightly exceeded and the warning can be ignored safely. You can print the orbital occupations for the orbitals affected by the +U correction with this print key at the PRINT_LEVEL medium to check the actual occupation. Note, that the U values found appropriate with PW codes in the literature are not necessarily optimal for CP2K, too. CP2K often gives a similar effect, e.g. impact on the band gap, for smaller U values.

 

HTH

 

Matthias

 

--
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/f41c9258-21be-4389-a69f-55aab2a808d4n%40googlegroups.com.

Nicholas Winner

unread,
Apr 21, 2021, 1:15:08 PM4/21/21
to cp2k
Thanks for the reply Matthias. I've had moderate success following what you suggested.

One issue that still persists, however, is that for some calculations I can only get convergence by using PLUS_U_METHOD LOWDIN. For example, MnO with Ueff=1.75eV, diagonalization and small smearing (T=100). I've tried using diagonalization and OT with MULLIKEN to get convergence, but sometimes LOWDIN seems to be the only way. Unfortunately, Lowdin does not support forces, so relaxations are out of the question. 

Do you have a suggestion for how to work around this?

-Nick

Krack Matthias (PSI)

unread,
Apr 21, 2021, 1:30:13 PM4/21/21
to cp2k

Hello Nick

 

The selection of the initial guess can help to achieve convergence. Could you provide a case which fails to converge with MULLIKEN? Otherwise, it is difficult to give further hints.

 

Matthias

 

Nicholas Winner

unread,
Apr 21, 2021, 1:36:02 PM4/21/21
to cp2k
Here is my input file for Mn3O4. 
cp2k.inp

Krack Matthias (PSI)

unread,
Apr 21, 2021, 2:12:25 PM4/21/21
to cp...@googlegroups.com

At a first glance, the basis set TZV2PX-MOLOPT-GTH-q6 for O is not well suited for condensed phase systems. You should rather use MOLOPT basis sets with “SR” in the name.

 

HTH

 

Matthias

 

Nicholas Winner

unread,
Apr 21, 2021, 6:58:19 PM4/21/21
to cp2k
So something very strange happened.

I'm hesitant to switch O to a SR basis, because only SZV/DZVP are available. SZV isn't production quality and DZVP really isn't near the basis set limit, but I gave it a try to see what would happen. The warning about negative DFT+U persists, but the energy changed to  -15383.653493175701442. I then re-initialized with TZV2PX and converged to  -15385.210534637011733. 

So to summarize I have three values for U=1.75
-15378.643602015281431 with TZV2PX
-15383.653493175701442 with DZVP
-15385.210534637011733 after restarting TZV2PX from the DZVP wfn.

These are enormous differences in energy and I don't know what to make of it. Any more ideas about what is going wrong based on my inputs?

Krack Matthias (PSI)

unread,
Apr 22, 2021, 4:04:33 AM4/22/21
to cp...@googlegroups.com

Hello Nick

 

The O DZVP-MOLOPT-SR basis set is accurate for condensed phase systems. If you want to explore the basis set limit for such systems, then you better use a plane waves code. You can add, of course, further polarization functions, but softer basis functions, as they are included in basis sets for molecular systems, will inevitably cause problems due to over-completeness which will result in ill-conditioned overlap matrices. The Mulliken orbital charges calculated with such overlap matrices are usually not reasonable and thus the calculation of the +U contribution is spoiled. The Löwdin analysis is more resilient than the Mulliken analysis in that respect, but it does not solve the basic problem. I am surprised that the restart using a TZV2PX basis set from a DZVP wfn restart file worked. That sounds more like garbage-in garbage-out.

 

Matthias

 

Matt W

unread,
Apr 22, 2021, 5:18:31 AM4/22/21
to cp2k
Just to note there are MOLOPT-TZVP-SR in the BASIS_MOLOPT_UCL files distributed with CP2K.
Matt

Nicholas Winner

unread,
Apr 22, 2021, 2:33:36 PM4/22/21
to cp2k

Matthias, after investigating further, it does seem like using the SR potential is fairly consistent in giving good results for the DFT+U calculations, even when Mulliken leads to a negative +U energy. I have previously not used DFT+U, and always used TZV2PX for O finding no problems, but the population analysis does seem to require a better conditioned overlap matrix. 

Matt mentioned TZV-SR basis sets in the UCL file, which is true, but not for O. In fact there are a few elements which only go up to DZVP-SR across all the BASIIS_MOLOPT_{, UCL, LnPP1} files. Some of these are O, C, Cl, F, H, etc. which are very common elements. H for instance only has SZV-SR and the rest are non-SR. If I believe that SR should always be used for stability in condensed phase, then it seems like these need some new basis sets, at least up to TZV quality. Maybe you disagree and DZVP is sufficient for these elements for some reason, while others really need larger bases?

Is there a reference for the procedure for fitting the SR basis sets? The molopt sets are described in "Joost VandeVondele and Juerg Hutter, J. Chem. Phys. 127, 114105 (2007)", but the BASIS_MOLOPT file simply says "variants of these basis sets using less and thus less diffuse primitives " were used for the SR versions, which is a little bit vague to me.

Krack Matthias (PSI)

unread,
Apr 22, 2021, 4:28:41 PM4/22/21
to cp...@googlegroups.com

Nick, we are talking about dense solid state systems. For such system, space is already well covered by a DZV basis set at each atom. The situation is, however, different already for molecular liquids like water or for systems with larger voids (MOFs, surfaces) or even small isolated molecular systems. Soft functions are required to describe the electron density decaying into the void regions. More polarisation function (second set of d and a set of functions) certainly improve the description, but  these function sets increase the computational costs significantly while their impact is rather moderate for DFT applied for O in dense solid state systems. I tried TZV/QZV and/or more polarisation functions for O in such systems, but the gain was rather small compared to the additional computational costs. There are tutorials showing how you may generate new or augment existing MOLOPT basis sets. It is also important to employ a balanced set of basis sets for a system which is, of course, also true for other system types. A mixing of SR and non-SR basis sets should always be done with care. The optimisation of the SR basis sets includes also the condition number of the overlap matrix as a parameter during the basis set generation procedure which allows for a control of its numerical stability.

Nicholas Winner

unread,
Apr 25, 2021, 9:12:22 PM4/25/21
to cp2k

Matthias, that all makes a lot of sense. I think I will do a few formation energy comparisons at the DZV / TZV levels for my own sanity, but that's for me to worry about.

The one thing that remains is the issue of Lowdin vs Mulliken. I still have some TMOs that struggle to converge even once I switch to the SR basis sets, but work when using Lowdin+SR basis sets.  Of course, Lowdin does not have forces implemented. Do you happen to know why they have yet to be implemented? Based on my knowledge, Lowdin and Mulliken are just two different ways to get q_{mu, nu}, and so I don't see why Lowdin would not be implemented and Mulliken is. If it is just a matter of effort, I might fork Cp2k and take a look at the source, but if there is a theoretical reason for it, then I won't bother.

Thanks for your help,
Nick

hut...@chem.uzh.ch

unread,
Apr 26, 2021, 3:42:16 AM4/26/21
to cp...@googlegroups.com
Hi

some information on MOLOPT basis sets and Lowdin charge derivatives:

Original MOLOPT basis sets were fully optimized for molecular energies
under a constraint to keep a low condition number. This optimization was
very time consuming and the resulting smallest exponents are leading
to very expensive collocate/integrate. Joost then generated a second set
(SR) with less primitives that keep the MOLOPT advantages but have
better CPU performance. Later, Florian Schiffmann introduced a wavefunction
fitting that reduced the optimization procedure.
I have recently generated new sets based on the SR model for PBE/SCAN/PBE0
pseudopotentials (see https://github.com/juerghutter/BASIS -> MOLOPT).
However, these are not yet tested.

Gradient for the Lowdin charges require the derivative of sqrt(S), whereas
Mulltiken only need derivatives of S. The task is to find a way to calculate
efficiently the derivatives of sqrt(S) as a function of the derivative of S.
For example, this can be done for S^-1 and is widely used in QC codes.
We have recently derived such a scheme for the gradient of the sTDA method.
So now what is missing is the implementation for "plus U".

best regards

Juerg Hutter

--------------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Institut für Chemie C FAX : ++41 44 635 6838
Universität Zürich E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: "cp2k" <cp...@googlegroups.com>
From: "Nicholas Winner"
Sent by: cp...@googlegroups.com
Date: 04/26/2021 03:12AM
Subject: Re: [CP2K:15237] DFT+U unphysical warning
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/df32ddfc-960a-4773-8dfa-82e9e4253bdan%40googlegroups.com.

Krack Matthias (PSI)

unread,
Apr 26, 2021, 7:22:52 AM4/26/21
to cp...@googlegroups.com
Nick, I have just committed a change to the current CP2K development (trunk) version which allows for the calculation of forces with the DFT+U method Lowdin and thus for GEO_OPT runs. The analytical stress tensor is still missing and thus CELL_OPT runs are still not possible. You may try it for your cases that did not work with Mulliken. Any feedback is welcome. We plan a CP2K 8.2 release in May which would include this change in case there are no problems.

While the computational overhead for the DFT+U method Mulliken is negligible, Lowdin requires the calculation of S**(1/2) for each new atomic configuration. The calculation of S**(1/2) requires currently a diagonalization of overlap matrix S, which can become a bottleneck for large systems. Moreover, S**(1/2) is a full matrix that does show the sparsity of S.

Matthias

> -----Ursprüngliche Nachricht-----
> Von: cp...@googlegroups.com <cp...@googlegroups.com> Im Auftrag von
> hut...@chem.uzh.ch
> Gesendet: Montag, 26. April 2021 09:42
> An: cp...@googlegroups.com
> Betreff: Re: [CP2K:15238] DFT+U unphysical warning
> https://groups.google.com/d/msgid/cp2k/66f4854a-2810-4a94-a5d7-
> 6372c2103c4cn%40googlegroups.com.
>
>
>
> --
> 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/a7c40bdf-fdfb-4448-a20e-
> bb24a88a7ee9n%40googlegroups.com.
>
> --
> 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/1f67b87b-7768-493b-b4e3-
> e5d46d426847n%40googlegroups.com.
>
> --
> 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/df32ddfc-960a-4773-8dfa-
> 82e9e4253bdan%40googlegroups.com.
>
>
> --
> 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/OFF29C6CBE.8C604B5F-
> ONC12586C3.0029E554-C12586C3.002A4F51%40lotus.uzh.ch.
Reply all
Reply to author
Forward
0 new messages