Variablecell does not converge

[mfherbst]

The attached cif file with this script:

using DFTK
using PseudoPotentialData
using GeometryOptimization
using AtomsIO

system = load_system("AlN_mp-661_conventional.cif")
pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
calc = DFTKCalculator(; model_kwargs = (;functionals=PBE(), pseudopotentials, symmetries=true),
                        basis_kwargs = (;  Ecut=30, kgrid=(4,4,4) ))

results = minimize_energy!(system, calc;
                           tol_forces=1.e-4, tol_virial=Inf, verbosity=2, variablecell=true, maxstep=Inf, )

does not converge. With variablecell=false it converges in 2 steps.

[mfherbst]

AlN_mp-661_conventional.txt

(renamed to txt because otherwise github complains)

[cortner]

before assuming this is a bug in the optimizer, can you

• (1) finite-difference check the virials for the DFT calculator?
• (2) After that, one should finite-difference check the mapping dofs -> energy. (for bugs in the (positions, cell) -> dofs mapping

Both of these can be tested in random configurations and in the configuration where the optimizer fails to make progress.

[WolframMiller]

using DFTK
using Plots
using PseudoPotentialData
using Unitful
using UnitfulAtomic

using AtomsIO        # Enables only Julia-based parsers

using PythonCall
using ASEconvert

using GeometryOptimization

using JLD2

using WriteVTK

function main()

    # Parameters

    E_cut = 400u"eV"

    k_grid = [4,4,4]

    temp = 1e-3
    
    pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
    
    tol_forces_A= 3.e-4  #  5.e-1u"eV/Å"
    tol_forces_B= 1.e-3
    tol_virial_B= 0.607

    #   geometry
    file_cif = "AlN_mp-661_conventional.cif"
    
    system = load_system(file_cif)
    println(system)

    num_nodes = length(system)  # number of nodes/particles
   

    #  Calculator for GeometryOptimization
    
    calc = DFTKCalculator(;
                          model_kwargs = (;functionals=PBE(), pseudopotentials, temperature=temp),
                          basis_kwargs = (;  Ecut=E_cut, kgrid=k_grid )
                          )
    ##

    ##-- Step 1

    # optimizers OptimLBFGS()  OptimCG() OptimCD()

    resultsA = minimize_energy!(system,calc,OptimCG(); tol_forces=tol_forces_A, verbosity=4, variablecell=false)
    
    system_stepA = resultsA.system  # relaxed system

    energyA = resultsA.energy
    system_stateA = resultsA.state

    model = model_DFT(system_stepA; functionals=PBE(), pseudopotentials=pseudopotentials, temperature=temp,
                      smearing=DFTK.Smearing.FermiDirac() )
    basis = PlaneWaveBasis(model; Ecut=E_cut, kgrid=k_grid )
    
    scfres = self_consistent_field(basis)
    scfres.energies

    
    f = compute_forces(scfres)
    println(f)

    
    # Save everything
    
    save_scfres("scfres_fixed.jld2",system_stateA)
    
    
    println(" -------  Step 1  finished   --------")

    unit_cell_unit = cell_vectors(system)
    unit_cell_unit_A = cell_vectors(system_stepA)
    
    unit_cell= Array{Float64,2}(undef,3,3)
    unit_cell_A= Array{Float64,2}(undef,3,3)

    lattice_constant   = Array{Float64}(undef,3)
    lattice_constant_A = Array{Float64}(undef,3)

    for i=1:3
        lattice_constant[i] = 0.0
        lattice_constant_A[i] = 0.0
        for j=1:3
            unit_cell[i,j] = ustrip(unit_cell_unit[i][j])
            lattice_constant[i] += unit_cell[i,j]*unit_cell[i,j]
            unit_cell_A[i,j] =  ustrip(unit_cell_unit_A[i][j])
            lattice_constant_A[i] += unit_cell[i,j]*unit_cell_A[i,j]
        end
        lattice_constant[i] = sqrt(lattice_constant[i])
        lattice_constant_A[i] = sqrt(lattice_constant_A[i])
    end

    alpha = 30/180*pi
    
    pos = Array{Float64}(undef,3)
    pos_A = Array{Float64}(undef,3)
    for l=1:4
        pos_with_units   = position(system,l)
        pos_without = ustrip(pos_with_units)
        pos[1] = (pos_without[1] +  pos_without[2]*tan(alpha) ) / lattice_constant[1]
        pos[2] = (pos_without[2] / cos(alpha)                 ) / lattice_constant[2]
        pos[3] = (pos_without[3]                              ) / lattice_constant[3]
        pos_with_units_A = position(system_stepA,l)
        pos_without_A = ustrip( pos_with_units_A)
        pos_A[1] = (pos_without_A[1] +  pos_without_A[2]*tan(alpha)) / lattice_constant_A[1]
        pos_A[2] = (pos_without_A[2] / cos(alpha)                ) / lattice_constant_A[2]
        pos_A[3] = (pos_without_A[3]                             ) / lattice_constant_A[3]
        println(" pos $l  coord  old  $(pos)   new $(pos_A)" )
    end
    
    println(lattice_constant_A)

    
    ##-- Step 2
    resultsB = minimize_energy!(system_stepA,calc,OptimCG(); tol_forces=tol_forces_B, tol_virial=tol_virial_B, verbosity=4, variablecell=true)
    println(" -------  Step 2  finished   --------")
    
    system_stepB = resultsB.system  # relaxed system
    system_stateB = resultsB.state
    
    unit_cell_unit_B = cell_vectors(system_stepB)
    unit_cell_B= Array{Float64,2}(undef,3,3)
    
    lattice_constant_B = Array{Float64}(undef,3)

    for i=1:3
        lattice_constant[i] = 0.0
        lattice_constant_B[i] = 0.0
        for j=1:3
            unit_cell[i,j] = ustrip(unit_cell_unit[i][j])
            lattice_constant[i] += unit_cell[i,j]*unit_cell[i,j]
            unit_cell_B[i,j] =  ustrip(unit_cell_unit_B[i][j])
            lattice_constant_B[i] += unit_cell[i,j]*unit_cell_B[i,j]
        end
        lattice_constant[i] = sqrt(lattice_constant[i])
        lattice_constant_B[i] = sqrt(lattice_constant_B[i])
    end



    
    pos_B = Array{Float64}(undef,3)

    for l=1:4
        pos_with_units   = position(system,l)
        pos_without = ustrip(pos_with_units)
        pos[1] = (pos_without[1] +  pos_without[2]*tan(alpha) ) / lattice_constant[1]
        pos[2] = (pos_without[2] / cos(alpha)                 ) / lattice_constant[2]
        pos[3] = (pos_without[3]                              ) / lattice_constant[3]
        pos_with_units_B = position(system_stepB,l)
        pos_without_B = ustrip( pos_with_units_B)
        pos_B[1] = (pos_without_B[1] +  pos_without_B[2]*tan(alpha)) / lattice_constant_B[1]
        pos_B[2] = (pos_without_B[2] / cos(alpha)                ) / lattice_constant_B[2]
        pos_B[3] = (pos_without_B[3]                             ) / lattice_constant_B[3]
        println(" pos $l  coord  old  $(pos)   new $(pos_B)" )
    end

    println(lattice_constant)
    println(lattice_constant_B)

    # Save everything
    
    save_scfres("scfres_fixed.jld2",system_stateB)
    
    return
    
end

• After some iterations no change in max(Force) and max(Virial)

• Then the error occur:
  ERROR: Mismatch in dimension between guess (772 and Hamiltonian (770)

• convergence is also slowed down to no progress for variablecell = false
  max(force) is about 2e-4

How can one check the
finite-difference check the mapping dofs -> energy ?

Could one print the the coordinates and lattice constants in every iteration step ?

[mfherbst]

How can one check the finite-difference check the mapping dofs -> energy ?

There are the functions eval_objective and eval_gradient in https://github.com/JuliaMolSim/GeometryOptimization.jl/blob/master/src/dof_management.jl, that's one would need to use.

Could one print the the coordinates and lattice constants in every iteration step ?

There is no simple way to do that (yet), but it can in principle be achieved using the callback parameter to minimize!. Perhaps taking a look at https://github.com/JuliaMolSim/GeometryOptimization.jl/blob/master/src/callbacks.jl helps to achieve that.

[Technici4n]

This is a classic problem for a plane-wave basis. Still, it took me a long time to figure it out.

The following code does converge:

using DFTK
using PseudoPotentialData
using GeometryOptimization
using AtomsIO

system = load_system("AlN_mp-661_conventional.cif")
pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
calc = DFTKCalculator(; model_kwargs = (; functionals=PBE(), pseudopotentials, symmetries=true,
                                          kinetic_blowup=BlowupCHV()),
                        basis_kwargs = (; Ecut=30, kgrid=(4,4,4) ))

results = minimize_energy!(system, calc;
                           tol_forces=1.e-4, tol_virial=Inf, verbosity=2, variablecell=true, maxstep=Inf, )

In particular note the kinetic_blowup=BlowupCHV() which is required to ensure that as the cell is resized the total energy remains continuous. Without it the number of basis functions can change abruptly such that the total energy is discontinuous, forcing the line search to pick ridiculously small steps (and causing other weird problems such as the gradients "jumping" due to being discontinuous)...

[cortner]

that sounds to me as if it is not an optimization problem but an issue with how plane wave DFT is implemented. I'm not a DFT expert so don't understand why you would resize the basis when the cell changes.

Naively I'd just add some kind of feedback/callback mechanism that terminates the optimization when the model changes so it can be restarted, then one could add an auto-restart wrapper.

But surely this has been solved 100 times over?

[Technici4n]

Indeed this is not an issue with this package, but rather with the plane-wave DFT package. I don't think that restarting the optimization is really practical because that does not solve the fundamental issue of the optimization landscape being bad. For your info here are some practical fixes from established codes:

1. Quantum ESPRESSO strategy: Fix the basis in the beginning of the optimization, and then optimize with a changing cell but fixed basis. Then once the optimization converged perform a final energy evaluation with the basis computed from the optimized cell. If there is a large energy difference, this suggests that the basis is too small anyway. Quite elegant IMO. (source: https://www.quantum-espresso.org/Doc/pw_user_guide/node11.html#SECTION00044030000000000000)
2. ABINIT strategy: Use a larger basis with a smooth penalty for the added basis functions. Also quite a neat approach IMO.

In DFTK strategy is 2 is implemented by adding a kinetic_blowup=BlowupCHV() in the right place. :)

[cortner]

ok thanks for the explanation.