Improved setup routines for charged systems

[umbriquse]

Looking through parts of the code, and documentation I did not notice a way to implement a given atom with a charge (e.g Fe^+2). If I missed how to do this in the docs or code I apologize!

I know it's possible to declare the number of electrons in the system, but in the case of using ions it didn't seem to know which atoms were missing the electrons.

I'm currently using a very rough way of implementing ions into a model, but I lack the required knowledge if this "fix" is appropriate.

function loadAtoms(elementSymbol::Union{Int,Symbol}; family::String = "hgh", ionCharge::Int = 0, core::Symbol = :fullcore, functional = "lda")::ElementPsp
    pspfile = list_psp(elementSymbol; family = family, functional = functional, core = core)
    pspfile = ifelse(length(psp) == 0, psp.identifier, psp[1].identifier)
    psp = load_psp(pspfile)

    @unpack Zion,rloc,cloc,rp,h,identifier,description = psp
    #TODO: This is very sketchy, I have no idea whether if changing the Zion should affect rloc,cloc...
    ElementPsp(elementSymbol; psp = PspHgh(Zion - ionCharge, rloc, cloc, rp, h; identifier = identifier, description = description))
end

If this is an inappropriate method for implementing ions and I should not trust the results, could you let me know?

[antoine-levitt]

Hi,

I know it's possible to declare the number of electrons in the system, but in the case of using ions it didn't seem to know which atoms were missing the electrons.

In the case of DFT charge is free to move around as it wants, so there's no specific notion of "this atom is missing two electrons", more "the total system is missing two electrons". So just changing the number of electrons should be what you need. Changing Zion is likely not a good idea, because you'll change the pseudopotential, which you probably don't want to do.

That said there are specific issues with isolated and charged systems which we haven't really considered yet, so be careful with that. What system do you want to do?

[mfherbst]

One thing that also came up in other projects was that setting the number of electrons explicitly instead of just saying something like "one extra electron" is a bit unhandy. I guess one could make a nice wrapper in the model_atoms function by just supporting an additional keyword argument charge or total_charge that takes the desired charge of the system.

Let me expand a bit on the warning of @antoine-levitt: Of course you can do calculations with charged unit cells in DFTK, but be careful that the periodic boundary conditions we always empose will cause artefacts from the charged defect centres of neighbouring unit cells repelling each other across the periodic boundary. For this reason you will need pretty large supercells to get a converged result. There are ways to correct for this and improve convergence wrt. supercell size (e.g. Makov-Payne corrections DOI 10.1103/physrevb.51.4014), but they are not yet available in DFTK.

[antoine-levitt]

Makov-Payne is for crystal defects, which is much more tricky. For isolated systems I don't really know what people do; possibly taking lots of vacuum is enough? We should look up what people do but it doesn't look like it's been so seriously tackled, there should be space for improvement on the state of the art here...

[umbriquse]

@antoine-levitt Thanks for the advice! The systems that I would like to simulate involve Fe-2 atoms interacting with metallic surfaces. Or, more generally, ionic atoms interacting with metallic surfaces. I know for abinit when calculating the surface energy of a material, they suggest to increase the vacuum until the effect the periodic lattice has on the system diminishes to an acceptable level. I imagine that concept can be generalized for isolated systems.

[mfherbst]

Yes, increasing the vacuum is amongst the things you have to do, but another problem is that the positively charged ions might repell each other across the domain you simulate (because we have no correction scheme implemented which would partly "cutoff" the Coulomb tail from one unit cell to the next). So you will also need to make sure that your metallic surface is sufficiently large. In practice this means that on top of increasing the vacuum size you also need to check convergence wrt. increasing the surface cell dimensions, which can get quite expensive. Using appropriate corrections you can make that convergence quite a bit faster, but as mentioned we don't have them implemented yet.

[antoine-levitt]

Cool, I'd be interested in seeing how that goes! In the case of a metallic surface hopefully screening is efficient enough that you don't need a huge surface. All the correction schemes I've seen are for semiconductors.

[mfherbst]

Same here :smile: Please let us know!

[stur86]

Something only tangentially connected to this, but is there a way to add a non-uniform charge background to the system, and/or change the charge and mass of the particle we solve for (replacing e.g. electrons with positrons or muons)? I'm asking because personally I'm interested in doing some two-component DFT and muonic systems, but while it's a pretty exotic application, it should also be quite easy to implement, in theory.

[mfherbst]

non-uniform charge background

You can add an arbitrary potential to the Hamiltonian, for example using the ExternalFromReal term (see this example for a way to add an electric field for example

change the charge and mass of the particle we solve

For the electrons the charge is implicitly used in the hamiltonian terms that govern how the electrons interact with the nuclei and how they act with each other. Of course if you adapt the terms appropriately (i.e. modify the local potentials generated by the nuclei and adapt the Hartree term) you can cheat the electron density to be a "muon density" or a "positron density".

(For the sake of keeping this issue focused on one topic could you please open a new issue if this does not answer your question).

[antoine-levitt]

Might get around to do it (for purely isolated systems), it's a good relaxing summer project.

1. Implement monopole correction for the electrostatics. Test system: H+, which hopefully should converge exponentially fast with system size
2. Implement dipole correction. Test system: H2O

API-wise we should have something in the model that tells us in which direction(s) the system is periodic and in which it is isolated. After this the question is what to do for partially periodic systems, but I think that's a research project.

[antoine-levitt]

Is H3O+ doable with DFT, maybe with collinear spin? Or are open shell molecules just not doable at that level?

[mfherbst]

Re dipole corrections. I have a draft branch with this I can share. It has a reasonable API. I think it should also work but it might have some strange unit issues.

[antoine-levitt]

Oh really? How do you do it? I was planning on subtracting a reference (gaussian) density with an analytic potential, like in the anyon stuff

[mfherbst]

No I wanted to try that, but never got to it. I follow GPAW. just compute the dipole moment in the cell (taking the centre as a reference) and from that compute a compensating field across the cell to get nice flat work functions.

[antoine-levitt]

Do you have a reference for that?

[mfherbst]

I think it's this one: http://dx.doi.org/10.1103/physrevb.59.12301

[antoine-levitt]

That's for surfaces, not molecules, right?

[mfherbst]

absolutely. Ah yes sorry you thought about molecules. I only read dipole corrections and heard my cue ;).

[antoine-levitt]

OK yeah that's more involved than what I was thinking about - let's talk about that later!

[mfherbst]

Sure. I might still pick it up again though ... not sure yet.

[antoine-levitt]

I had a look at the paper, it looks like a reasonable thing to do, but it's not quite what I had in mind - in the anyons stuff I use a reference gaussian density of finite width having the same characteristics (monopole/dipole) as the actual density. Here they essentially take only the dipole but don't have a finite width. I need to think about it some more and play with it in molecules.

[antoine-levitt]

OK, I had a closer look at the paper, I think adding and subtracting a gaussian is better. I'll code it up (should be a small modification of the scheme for molecules), it would be great to compare to this approach if you've got it coded up

[antoine-levitt]

OK thinking about it more this is non-trivial because it also couples with how the local potential from nuclei is constructed. We should think about it more seriously, meaning it'll have to wait a bit.