Spin unrestricted calculations with DFTpy

In this tutorial we will perform spin unrestricted calculations with OF-DFT using DFTpy

Note: Install the dev branch of DFTpy.

import numpy as np
from ase.build import bulk
import matplotlib.pyplot as plt

First we need to load the necessary modules form DFTpy

from dftpy.api.api4ase import DFTpyCalculator
from dftpy.config import DefaultOption, OptionFormat
from dftpy.functional import LocalPseudo, Functional, TotalFunctional
from dftpy.optimization import Optimization
from dftpy.ions import Ions
from dftpy.field import DirectField
from dftpy.grid import DirectGrid
from dftpy.math_utils import ecut2nr
from dftpy.formats import io

def scale_density(rho, m):
    if rho.rank != 2:
        raise Exception("Rho must be rank 2")
    nelec = rho.integral()
    nnelec = nelec + np.array([m/2.0,-m/2.0])
    rho[0] *=  nnelec[0]/nelec[0]
    rho[1] *=  nnelec[1]/nelec[1]
    return rho

Build the FCC Al crystal structure with ASE and then change the units to DFTpy units

atoms = bulk('Al', 'fcc', a=4.05)
ions = Ions.from_ase(atoms)

Define the Exchange-correlation, Kinetic energy and Hartree functionals

XC = Functional(type='XC',name='LDA', libxc=False)
HARTREE = Functional(type='HARTREE')
KEDF = Functional(type='KEDF', name='TFvW', y=1)
opt_options = {'econv' : 1e-7*ions.nat, 'maxiter': 50}

Build the grid and define the Local pseudopotential

path = "../DATA"
PP_list = {'Al':'al.lda.recpot'}

nr = ecut2nr(ecut=90, lattice=ions.cell)
grid = DirectGrid(lattice=ions.cell, nr=nr)
PSEUDO = LocalPseudo(grid=grid, ions=ions, PP_list=PP_list)
rho = DirectField(grid=grid,rank=2)
rho[:] = ions.get_ncharges() / ions.cell.volume / rho.rank

setting key: Al -> al.lda.recpot

Changing the number of beta and alpha lectrons

rho = scale_density(rho,0.1)

rho.integral()

array([1.55, 1.45])

Perform the optimization of the density

evaluator = TotalFunctional(KE=KEDF, XC=XC, HARTREE=HARTREE, PSEUDO=PSEUDO)
opt = Optimization(EnergyEvaluator=evaluator, optimization_method='TN', optimization_options=opt_options)

rho = opt.optimize_rho(guess_rho=rho)

Step    Energy(a.u.)            dE              dP              Nd      Nls     Time(s)
0       -2.048930037247E+00     -2.048930E+00   2.989137E-01    1       1       2.462769E-02
1       -2.107897495307E+00     -5.896746E-02   4.295104E-02    4       40      4.644361E-01
2       -2.111292714918E+00     -3.395220E-03   3.350295E-03    28      40      1.114891E+00
3       -2.111412636599E+00     -1.199217E-04   2.943865E-04    20      40      1.690181E+00
4       -2.111424538274E+00     -1.190167E-05   2.342158E-05    28      19      2.133016E+00
5       -2.111425031454E+00     -4.931807E-07   1.915958E-06    18      19      2.487958E+00
6       -2.111425095467E+00     -6.401236E-08   1.582082E-07    24      19      2.898481E+00
7       -2.111425100518E+00     -5.051136E-09   6.574188E-08    22      19      3.287766E+00
#### Density Optimization Converged ####
Chemical potential (a.u.): [0.28534709 0.28956009]
Chemical potential (eV)  : [7.76468977 7.87933132]

Finally compute the total energy of the system

evaluator.Energy(rho)

-2.1114251005178186