I am trying to apply Hubbard U (DFT+U) to Ce in cerium based oxides including CeO2 and Ce2O3. However, as you see from header of Ce psp (see [I] below), our Ce psp (RRKJ, generated by OPIUM) is designed to have zero 4f occupation. When I applied +U to Ce using QE V4.3.1, I got following error message:
from offset_atom_wfc : error # 9 wrong offset So, I searched previous threads and I figured out that I needed to modify QE code (http://www.democritos.it/pipermail/pw_forum/2010-July/017622.html). I did modification and recompiled QE code. It looked fine but I found weird energetics for reduction of CeO2 to Ce2O3 in comparison with data from PRB 75, 035109 (2007). So, I looked into output files and found out that occupation matrix for both CeO2 and Ce2O3 are somewhat unphysical (at least for me): *U_eff=0.01 eV Tr[ns(na)] CeO2 = 1.1690424 Ce2O3 = 1.3908759 *U_eff=3.0 eV Tr[ns(na)] CeO2 = 0.9963736 Ce2O3 = 1.2650403 *U_eff=9.0 eV Tr[ns(na)] CeO2 = 0.5972449 Ce2O3 = 1.1056369 I thought that even with non-orthogonal terms, at least Tr[ns(na)] for CeO2 should be very close to 0.0 as Ce exists as Ce4+ in CeO2 (with ionic picture) and should have empty 4f state. With analyzing Ce psp itself, I wanted to know whether the modification of QE discussed above is safe or not. Answer of Dr. Sclauzero to my question was as follows: "Briefly, I think that your problem stems exactly from that fix which I proposed and is now included in the latest PWscf version. That fix will actually introduce a problem if the atomic wave function which you want to use for making projections have zero occupations. I thought that this would never be the case, but probably one should not exclude it a priori. A further problem arises if you generate such a pseudopotential (Ce with zero 4f occupation) and want to use it for making the LDA+U projections on the empty states, since the LD1 code usually saves in the PP file only the wavefunctions with nonzero occupations. First, you should check if you PP has the 4f wavefunctions. If it does, you can remove the check upf(nt)%oc(n) > 0.D0 as a quick fix. If not, then you need to regenerate the PP". Based on his comment, I did some test calculations after modifying our Ce psp--simply made zero 4f occupation as nonzero, 0.01 (In (I), 430 0.00 --> 430 0.01). However, there is no total E change between the calculation with Ce psp of zero 4f occupation and that with nonzero 4f occupation (E difference <0.001/CeO2). Also, change in Tr[ns(na)] is minute (e.g. 0.9938046 --> 0.9963736) My questions to QE users are: - Maybe I need to fill 4f more? - Do you have any explanations about why is occupation matrix sum close to 1.0 not 0.0? Again, while there are some non-orthogonal terms, I think that for CeO2 Tr[ns(na)] for Ce should be close to zero not one. Did I miss something? - How about using starting_ns_eigenvalue? By the way, while I tried to use this, I couldn't understand instruction. For instance, in example below ([II]), if I want to make initial occupations as zero, how can I do? In starting_ns_eigenvalue(m,ispin,I), I can see that I should have starting_ns_eigenvalue(m,1,1) and starting_ns_eigenvalue(m,2,1) if U is applied to spin-up and spin-down for atomic species 1. Then how about m? m=1,2,3,4,5,6, and 7? If not, what will be desirable values for m? Thank you for your advice and comments in advance. P.S. I post this while there is existing post regarding the questions above because: 1) The title of the existing post is not suitable and will never be easy for QE users to search. 2) I am not comfortable with adding threads under the post which contains my personal emails (not for the forum) which violate my privacy. ++++++++++++++++++++++++++++++ [I] [Atom] Ce 15 # norb: number of orbitals 100 2.00 -643.0 # nlm occ eigen(- means auto-generate) 200 2.00 -77.4 210 6.00 -68.4 300 2.00 -10.0 310 6.00 -6.0 320 10.0 -6.4 400 2.00 -6.4 410 6.00 -6.0 420 10.0 -6.0 500 2.00 -2.4 510 6.00 -2.0 600 0.00 -1.0 610 0.00 -1.0 520 0.00 -1.0 430 0.00 -1.9 +++++++++++++++++++++++++++++++ [II] U( 1) = 3.0000 U( 2) = 0.0000 alpha( 1) = 0.0000 alpha( 2) = 0.0000 atom 1 Tr[ns(na)]= 2.0000000 atom 1 spin 1 eigenvalues: 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 eigenvectors 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 5 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 7 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 occupations 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.143 ------- Chan-Woo Lee, Ph.D. Department of Chemistry University of Pennsylvania 231 South 34th Street Philadelphia, PA 19104-6323 Phone: 1-215-898-3564 (Office) Email: <mailto:leechanw at sas.upenn.edu> leechanw at sas.upenn.edu / <mailto:cwandtj at gmail.com> cwandtj at gmail.com -------------- next part -------------- An HTML attachment was scrubbed... URL: http://www.democritos.it/pipermail/pw_forum/attachments/20111214/f3815e8a/attachment-0001.htm
