DFTK.jl
DFTK.jl copied to clipboard
JuliaMolSim
WIP: Troullier-Martins pseudopotential parser
So this is my progress as of now. Sorry it has taken this long, I had no idea how pseudopotentials worked.
The parser is fully working, I am now stuck on where to go from here.
Looking at #397 @antoine-levitt suggested to
compute the integral of a radial function specified on a grid against a bessel function j(q r), for several q
But this format provides radial values for the pseudopotentials, and the first and second projector values. So I'm a bit confused as to which to use and also what this function is intended to do.
EDIT: @antoine-levitt would you perhaps be referring to equation 19 in the "Efficient pseudopotentials for plane-wave calculations" paper?
Great, thanks for this and happy new year! Yes indeed 15 and 19 are what I was thinking about. The psp file gives you the radial part on a grid of the local and nonlocal potentials (or possibly the KB form of the nonlocal pot directly, I'm not sure), and DFTK uses the matrix elements in a plane wave basis (psphgh.jl has both real and recip space for completeness, but only the recp space expressions are actually used in practice - at least for now). These matrix elements involve Fourier transforms of radial times spherical harmonics, which in turn are expressed as integrals of the radial part against bessel functions (see https://juliamolsim.github.io/DFTK.jl/dev/advanced/useful_formulas/). There are some scattered docs in the psphgh and in local.jl and nonlocal.jl (in terms/), if anything is unclear please let me know and I can explain or add docs.
To be clearer : since you're talking about projectors it means it's directly in KB form, the link with the TM paper is that the projectors are phiPP * Vnonloc (modulo normalisations and trickiness like that).
After looking at the code in nonlocal.jl it appears we have been pretty careful about separating generic pseudopotentialness and specific Hgh (@mfherbst to thank for that!), so essentially you just have to use a new type instead of psphgh and add the missing methods that pop up.
Hi @umbriquse and happy new year!
Thanks for giving this a go! Getting the HGHs to work was tedious enough so we never really bothered going beyond it, but TM pseudos are pretty wide-spread so definitely good to have them.
Indeed to get a new pseudo to work in DFTK you only need to slot it into the psp variable of the ElementPsp struct and make sure that all the required functions are implemented for PspTM. To be precise that is
eval_psp_projection_radial(KB Projector in Fourier space, i.e. equation 19)eval_psp_energy_correction(Correction of the ewald energy due to the PP ... have not found the expression on first sight in TM paper, but cannot be hard to compute it)eval_psp_local_real(local part in real space, currently unused, but should still be implemented for later)eval_psp_local_fourier(local part in fourier space the Vlocal(G - G') perhaps modulo different normalisation)
Indeed it's as simple as that. We also use psp.h and psp.lmax but that should be it.
I'll make a PR adding an abstract supertype and outlining what needs to be implemented, hopefully that'll clarify things.
Good idea. When researching this it occurred to me as well that this could be useful.
Hm, should we merge ElementPsp and PspHgh (ie have PspHgh inherit Element directly)?
I'm not sure because you might want to have an element of one sort but combine it deliberately with the pseudo of another element. If element and pseudo are kind of the same thing that might become confusing (e.g. what should the parser return? It doesn't know the Z only the Zion. How do you change the Z / element symbol after the parser has given you an Element object (currently it is not mutable)?)
The extra indirection is a bit confusing so it'd be good to get rid of it.
Because you want to play with the number of valence electrons and the charge?
Sure give it a try. But what I like about the current design is that you can just say ElementPsp(psp=load_psp(identifier)) and this does the right thing for every psp identifier (once it's implemented). So you don't need to specify the Psp type twice (once in the identifier and once in the type of element you use). That would be my main point of concern for the design if you merge both things.
Because you want to play with the number of valence electrons and the charge?
I'm not sure it makes sense. Let's discuss it sometimes
Sure give it a try. But what I like about the current design is that you can just say ElementPsp(psp=load_psp(identifier)) and this does the right thing for every psp identifier (once it's implemented). So you don't need to specify the Psp type twice (once in the identifier and once in the type of element you use). That would be my main point of concern for the design if you merge both things.
That's surface API, it shouldn't interfere with the design of the structs. It can be done by an external function (or constructor), no? I'm not going to do it for now, but something to keep in mind.
It can be done by an external function (or constructor), no?
I'm sure it can be done somehow. It's just something to keep in mind because types are to some extent more static than strings ... and the list_psp / load_psp stuff pops up at quite some places in user code, so we should invest some thought how changes can be least painful. Let's see when you have a design.
Yeah I'm not too confortable around the load/listpsp stuff, so I'll let you handle it if you're OK with that
I think the comment p(q) = ∫_{R+} r^2 p(r) j_l(q r) dr in eval_psp_projection_radial is missing a 4\pi (like for the local potential). Do you agree?
@umbriquse OK so with the new API in https://github.com/JuliaMolSim/DFTK.jl/pull/401 you just need to create a new subtype of NormConservingPsp and implement the documented methods.
I'm going to need someone to check over the 'eval_psp_projection_radial' function. I'm not confident that I've got it 100% correct.
I think you're still working off the old code without the refactoring PR, eg eval_psp_projection_radial has changed name to eval_psp_projector_fourier I think.
From a very cursory glance it looks like it should be doing the right thing. I think a very good starting point to compare is the @testset "Numerical integration to obtain fourier-space projectors" begin in tests/PspHgh.jl (which should also be extended to the pair (eval_psp_local_real, eval_psp_local_fourier)). After that, once you're confident the _fourier and _real are consistent, I would just try to set up a silicon computation, compare against another DFT and hope for the best!
compare against another DFT and hope for the best!
I think ABINIT has TM pseudos, which has the advantage that we have parsers for their densities (i.e. you could converge the SCF in ABINIT and really check the resulting SCF state term by term) ... but I would not do this detailed comparison unless you really need it for debugging, but good to have it in the back as an option.
@antoine-levitt From memory ABINIT has since discontinued the use of TM pseudos due to their inaccuracy.
EDIT: Here's the page it's under 'Norm-conserving pseudopotential tables'
Well the nice thing is that the functionality should really be the same for TM and ONCVPSP so you might as well just read the ONCVPSP pseudos?