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.