How to set up models

[killah-t-cell]

First, sorry for asking really basic questions about DFTK and DFT. I am super thankful for any help, and I'd understand if you'd prefer to just close the issue and let me figure it out on my own.

I am trying to set up a Pt system in DFTK, yet I can't seem to do so with convergence just yet (as in, I get SCF not converged. after 100 iters). I am often getting errors such as Unable to find non-fractional occupations and Ortho(X,Y) is failing badly. I am not sure whether this has to do with how I set positions (which I am doing fairly stochastically), kgrid, Ecut, or something else. But I wanted to share my code and see if you know.

# Si-tutorial example
a = 5.431u"angstrom"        
lattice_ = a / 2 * [[0 1 1.];
                   [1 0 1.];  
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms     = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
model = model_LDA(lattice_, atoms, positions)
kgrid = [4, 4, 4]     # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 7              # kinetic energy cutoff
basis = PlaneWaveBasis(model; Ecut, kgrid)
scfres = self_consistent_field(basis, tol=1e-8);
# this works!

# Pt–my code
a = 3.8u"angstrom"  
lattice_fcc = lattice(:fcc, a)
Pt = ElementPsp(:Pt, psp=load_psp("hgh/lda/Pt-q10"));
atoms = [Pt]
positions = [[1/2, 1/2, 1/2]]
model = model_LDA(lattice_fcc, atoms, positions)

# try 1
basis = PlaneWaveBasis(model; Ecut=15)
scfres = self_consistent_field(basis, tol=1e-8);
# LoadError: Unable to find non-fractional occupations that have the correct number of electrons. You should add a temperature.
# is kgrid always [1,1,1] for single atoms?

# try 2
basis_with_k = PlaneWaveBasis(model; Ecut=10, kgrid=[1,1,1])
scfres_with_k = self_consistent_field(basis_with_k, tol=1e-8);
# Warning: SCF not converged.

# try 3
basis_low_cut = PlaneWaveBasis(model; Ecut=3, kgrid=[1,1,1])
scfres_low_cut = self_consistent_field(basis_low_cut, tol=1e-8);
# LoadError: Ortho(X,Y) is failing badly, this should never happen

# Al-system following tutorial https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/lattice_constants/lattice_constants.html
a = 4.04
lattice_fcc = lattice(:fcc, a)
Al = ElementPsp(:Al, psp=load_psp("hgh/lda/Al-q3"));
atoms = [Al]
positions = [ones(3)]
model = model_LDA(lattice_fcc, atoms, positions)
basis = PlaneWaveBasis(model; Ecut=200, kgrid=[8,8,8])
scfres = self_consistent_field(basis, tol=1e-8);
# LoadError: 3 electron cannot be attained by filling states with occupation 2. Typically this indicates that you need to put a temperature or switch to a calculation with collinear spin polarization.

My ultimate goal is to "perform calculations to determine whether Pt prefers the simple cubic, fcc, or hcp crystal structure", so I am trying to write code to minimize the energy with varying a. The optimization part seems to work, but I unfortunately cannot get the Pt system to converge at all, so I am not getting correct results.

kgrid = [1, 1, 1]       # k-point grid
Ecut = 10                # kinetic energy cutoff in Hartree
tol = 1e-8              # tolerance for the optimization routine
Pt = ElementPsp(:Pt, psp=load_psp("hgh/lda/Pt-q10"));
atoms = [Pt]

function compute_scfres(a; lat=:sc)
    model = model_LDA(lattice(lat, a), atoms, [ones(3)])
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    scfres = self_consistent_field(basis; tol=tol / 10)
    return scfres.energies.total
end;

a_res = optimize(compute_scfres, 0.0, 7.0) 

So in conclusion, I'd love to hear some more about:

• How does one set positions for a optimization problems and generally?
• How can one wisely vary kgrid/set it correctly for single and multiatom cells?

Thank you!

Ps: this could be a really great tutorial to add to the docs once it is done. GPAW has a similar doc

[antoine-levitt]

First, sorry for asking really basic questions about DFTK and DFT. I am super thankful for any help, and I'd understand if you'd prefer to just close the issue and let me figure it out on my own.

No problem, we're super keen to have user feedback on what's clear and what's not!

Unable to find non-fractional occupations

That tells you to add temperature, and you should do that. Some codes have temperature by default, we prefer to be on the explicit side and don't have it, but it's mandatory to have temperature (or interpolation schemes) for metals

Ortho(X,Y) is failing badly, this should never happen

OK, that error message may not be the most helpful :-) I've seen this happen when 1) NaN are produced or 2) the plane-wave cutoff is very small. I'll try to make the error message better here. I can't reproduce your example (because it doesn't know lattice(:fcc, a)), would you mind adding println("$N $M") to the beginning of the LOBPCG function in src/eigen/lobpcg_hyper_impl.jl and posting the result here? (if you've never edited a package before, do dev DFTK in Pkg mode, edit the file in ~/.julia/dev/, and restart Julia (or use Revise))

is kgrid always [1,1,1] for single atoms?

If you want to model a single atom, yes. But in that case you should have a large box to avoid boundary effects. Kgrid is there to perform (implicitly) a supercell computation (repeat the unit cell kgrid times), it's necessary for bulk (otherwise you're modelling a single cell with periodic boundary conditions, which is not physically correct) but doesn't really help for isolated systems.

How does one set positions for a optimization problems and generally?

For highly symmetric structures (fcc/bcc/etc) the symmetry is fixed and there is only a single parameter a. So you pick a reasonable initial value of a and let the optimizer do its thing. For non-symmetric structures, or if you don't want to impose any symmetry assumptions, you can pick a (reasonable...) initial lattice and let it relax with a multivariable optimizer (using stresses to compute the derivatives)

How can one wisely vary kgrid/set it correctly for single and multiatom cells?

This is the nasty part about DFT that ideally we'd like to automate. In general you have to do a convergence study: start with a reasonable grid, pick a property (eg energy), increase the grid until the property doesn't change too much and use that value. You have to do this study (separately) for kgrid, Ecut (and, in theory, temperature, though in practice you can often fix a reasonably small value and not care too much about it) Eg abinit has a nice tutorial on how to perform a convergence study.

Ps: this could be a really great tutorial to add to the docs once it is done. GPAW has a similar doc

GPAW has absolutely stellar docs. We should do better on docs, but that takes time... Note that (for now at least) DFTK is intended to be more developer-friendly than end-user-friendly, which explains some of the design choices and docs, but we're very open to making the docs friendlier.

[killah-t-cell]

@antoine-levitt here's the lattice() function I wrote so you can reproduce the example

function lattice(type, a) 
    if type === :fcc
        return a / 2 * [[1. 1. 0]; [0 1. 1.]; [1. 0 1.]]
    end

    if type === :bcc
        return a / 2 * [[1. 1. 0]; [0 1. 1.]; [1. 0 1.]]
    end

    if type === :sc
        return a * I(3)
    end
end

[killah-t-cell]

Ok, adding a temperature did it, and the system converges now. Is using temperature around 0.01 usually what is advised? (I see you mentioned a "reasonably small value") and then I copied this from the metallic systems example.

kgrid = [1, 1, 1]       # k-point grid
Ecut = 300u"eV"               # kinetic energy cutoff in Hartree
tol = 1e-8              # tolerance for the optimization routine
temperature = 0.01
Pt = ElementPsp(:Pt, psp=load_psp("hgh/lda/Pt-q10"));
atoms = [Pt]

function compute_scfres(a; lat=:sc)
    model = model_LDA(lattice(lat, a), atoms, [zeros(3)]; temperature)
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    scfres = self_consistent_field(basis; tol=tol / 10)
    return scfres.energies.total
end;

a_res = optimize(compute_scfres, 3.0, 5.0) 
a0 = Optim.minimizer(a_res)
auconvert(Å, a0)

would you mind adding println("$N $M") to the beginning of the LOBPCG function in src/eigen/lobpcg_hyper_impl.jl and posting the result here?

The output was 5 12. I hit this error sometimes depending on what value I use for a or Ecut

kgrid = [1, 1, 1]       # k-point grid
Ecut = 300u"eV"               # kinetic energy cutoff in Hartree
tol = 1e-8              # tolerance for the optimization routine
temperature = 0.01
Pt = ElementPsp(:Pt, psp=load_psp("hgh/lda/Pt-q10"));
atoms = [Pt]

function compute_scfres(a; lat=:fcc)
    model = model_LDA(lattice(lat, a), atoms, [zeros(3)]; temperature)
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    scfres = self_consistent_field(basis; tol=tol / 10)
    return scfres.energies.total
end;

compute_scfres(2) # LoadError: The size of the plane wave basis is 1, and you are asking for 12 eigenvalues. Increase Ecut.
compute_scfres(4) #Ortho(X,Y) is failing badly, this should never happen
compute_scfres(5) # -24.69861613674 (converges in 9 steps)

DFTK is intended to be more developer-friendly than end-user-friendly

For what is worth, I am having a lot of fun with it. I like that things are made explicit, and is given me a lot of opportunities to learn about DFT.

[antoine-levitt]

The output was 5 12

Then it should have printed a warning, didn't it?

[killah-t-cell]

@antoine-levitt yes, I got a

┌ Warning: Your problem is too small, and LOBPCG might
│         fail; use a full diagonalization instead

[antoine-levitt]

OK, good. We should probably just turn that into an error.

[antoine-levitt]

To clarify: this means that your discretized Hamiltonian is of size 5 (ie the eigenfunctions are expanded on 5 plane waves), and DFTK needs 12 eigenstates to do its job, which obviously doesn't work.

[killah-t-cell]

Yeah, that makes sense! This was super helpful. Hopefully I can take it from here :) thanks a lot, Antoine

[antoine-levitt]

https://github.com/JuliaMolSim/DFTK.jl/pull/667

[mfherbst]

Ok, adding a temperature did it, and the system converges now. Is using temperature around 0.01 usually what is advised? (I see you mentioned a "reasonably small value") and then I copied this from the metallic systems example.

Just to answer this one: 0.01 is a bit on the large side. Sometimes this is ok, but a safer default is 1e-3.

[mfherbst]

Closing that for now. Feel free to reopen @killah-t-cell if you have more questions.