Ctrl+K

Density Optimization#

Introduction#

The core functionality of any orbital-free density functional theory code is to be able to perform density optimizations. In theory, we seek to minimize the energy functional subject to the constraints that the electron density is positive and normalized. We can hence form a Lagrangian,

\[\mathcal{L}[n] = E[n] - \mu \left[\int n(\mathbf{r}) d^3\mathbf{r} - N_e \right],\]

which yields an Euler-Lagrange equation

\[\mu - \frac{\delta E}{\delta n(\mathbf{r})} \Bigg\vert_{n_\text{gs}} = 0,\]

to be solved for the ground state density. While introduced as a Lagrange multiplier constant, \(\mu\) can be interpreted as the chemical potential based on the Euler-Lagrange equation.

PROFESS-AD adopts a direct minimization scheme to solve the Euler-Lagrange equation. To account for the density positivity and normalization constraints during the optimization, PROFESS-AD minimizes the energy with respect to an unconstrained variable \(\chi(\mathbf{r})\) instead, which is related to the constrained density \(n(\mathbf{r})\) via

\[n(\mathbf{r}) = \frac{N_e~\chi^2(\mathbf{r})}{\int~d^3\mathbf{r}'\chi^2(\mathbf{r}')},\]

where \(N_e\) is the number of electrons. The convergence condition is hence

\[\frac{\delta E}{\delta \chi} \Bigg\vert_{\chi_\text{gs}} = 0\]

By default, PROFESS-AD utilizes a PyTorch-based LBFGS optimizer to perform density optimizations. The default convergence criteria is that the energy difference after each iteration must be below 1e-7 eV (controlled by ntol), 3 times in a row (controlled by n_conv_cond_count). Other parameter options can be found in the description of the system.optimize_density() method in the System page.

Basic Example#

Let us consider a basic example with the following code

# set energy terms and functionals to be used
terms = [IonIon, IonElectron, Hartree, WangTeter, PerdewBurkeErnzerhof]

# use "get cell" to get the lattice vectors and fractional ionic coordinates of
# a face-centred cubic (fcc) lattice
box_vecs, frac_ion_coords = get_cell('fcc', vol_per_atom=24.8, coord_type='fractional')

# defining the ions in the system
ions = [['Al', 'al.gga.recpot', frac_ion_coords]]

# set plane-wave cutoff at 2000 eV
shape = System.ecut2shape(2000, box_vecs)

# create an fcc-aluminium system object
system = System(box_vecs, shape, ions, terms, units='a', coord_type='fractional')

# perform density optimization (by default n_verbose is False, but we want to
# display the progress of the density optimization)
system.optimize_density(ntol=1e-7, conv_target='dE', n_method='LBFGS', n_verbose=True)

# check the measures of convergence
dEdchi_max = system.check_density_convergence('dEdchi')
mu_minus_dEdn_max = system.check_density_convergence('euler')

print('Convergence check:')
print('Max |𝛿E/𝛿χ| = {:.4g}'.format(dEdchi_max))
print('Max |µ-𝛿E/𝛿n| = {:.4g}'.format(mu_minus_dEdn_max))

This sets up a face-centred cubic (fcc) aluminium system and performs a density optimization. With the default settings, termination of the density optimization is based on the energy differences dE. We can also use the other measures of convergence (which correspond to the extent that the Euler-Lagrange equations are obeyed) to make sure that the density has been optimized to the desired extent.

Starting density optimization
  Iter      E [eV]      dE [eV]       Max |𝛿E/𝛿χ|       Max |µ-𝛿E/𝛿n|
   0       9.677157        0            0.865837           3.23348
   1       8.403218     -1.27394        0.473481           2.20042
   2       7.702040    -0.701178        0.550324           3.37185
   3       7.548173    -0.153867        0.211747           1.38436
   4       7.507185    -0.0409875       0.107438           0.771551
   5       7.485764    -0.0214207       0.089752           0.769107
   6       7.458289    -0.0274754      0.0758372           0.628908
   7       7.429842    -0.0284467      0.0384012           0.223996
   8       7.421262   -0.00857984      0.0292648           0.124719
   9       7.418701   -0.00256096      0.0197278           0.106185
   10      7.417557   -0.00114492      0.0106644          0.0758516
   11      7.417031   -0.000525209     0.00736959         0.0539031
   12      7.416772   -0.000259712     0.00520798         0.0398683
   13      7.416638   -0.000133214     0.00248938         0.0136218
   14      7.416570   -6.88236e-05     0.00358424          0.018048
   15      7.416518   -5.1554e-05      0.00181108         0.00946193
   16      7.416493   -2.47783e-05     0.00164527         0.00834793
   17      7.416479   -1.46257e-05     0.00139808         0.00779606
   18      7.416465   -1.38169e-05     0.00105684         0.0047043
   19      7.416458   -6.52723e-06     0.00103864         0.00740266
   20      7.416455   -3.72415e-06    0.000876987         0.00351695
   21      7.416452   -2.1472e-06     0.000503422         0.00219884
   22      7.416451   -9.56176e-07    0.000351898         0.00153222
   23      7.416451   -2.24164e-07    0.000156075        0.000990379
   24      7.416451   -2.25588e-08    0.000131899        0.000806672
   25      7.416451   -1.69511e-08    0.000139933        0.000889683
   26      7.416451   -1.92808e-08    9.49472e-05        0.000544465
Density optimization successfully converged in 26 step(s)

Convergence check:
Max |𝛿E/𝛿χ| = 9.424e-05
Max |µ-𝛿E/𝛿n| = 0.0005445

Example with Custom Potentials (Quantum Harmonic Oscillator)#

Let us consider a use case where we want to use a custom external potential for the density optimization. The procedure is similar to the above, just that we have to supply a dummy ions parameter to initialize the system object before setting the potential to our desired one, remembering to change the electron number of the system too.

# the single electron quantum harmonic oscillator (QHO) is a non-interacting
# single-orbital system - hence, it can be modelled well with just the
# ion-electron interaction and Weizsaecker terms
terms = [IonElectron, Weizsaecker]

# use a large box to simulate such localied systems with periodic
# boundary conditions so that the electron density will approach zero
# at the box boundaries
L = 20.0
box_vecs = L * torch.eye(3, dtype=torch.double)

# set low energy cutoff of 300 eV
shape = System.ecut2shape(300, box_vecs)

# as we will set the external potential ourselves later, we just need to
# submit a dummy "ions" parameter (the recpot file and ionic coordinates
# are arbitrary for this example)
ions = [['-', 'al.gga.recpot', torch.tensor([[0.5, 0.5, 0.5]]).double()]]

# create system object
system = System(box_vecs, shape, ions, terms, units='b', coord_type='fractional')

# as we have used an arbitrary recpot file, we need to set the electron number explicitly
system.set_electron_number(1)

# QHO quadratic potential
k = 10
xf, yf, zf = np.meshgrid(np.arange(shape[0]) / shape[0], np.arange(shape[1]) / shape[1],
                         np.arange(shape[2]) / shape[2], indexing='ij')
x = box_vecs[0, 0] * xf + box_vecs[1, 0] * yf + box_vecs[2, 0] * zf
y = box_vecs[0, 1] * xf + box_vecs[1, 1] * yf + box_vecs[2, 1] * zf
z = box_vecs[0, 2] * xf + box_vecs[1, 2] * yf + box_vecs[2, 2] * zf
r = np.sqrt(x * x + y * y + z * z)
qho_pot = 0.5 * k * ((x - L / 2).pow(2) + (y - L / 2).pow(2) + (z - L / 2).pow(2))

# set external potential to QHO potential
system.set_potential(torch.as_tensor(qho_pot).double())

# perform density optimization
system.optimize_density(ntol=1e-7, n_verbose=True)

# compare optimized energy and the ones expected from elementary quantum mechanics
print('Optimized energy = {:.8f} Ha'.format(system.energy('Ha')))
print('Expected energy = {:.8f} Ha'.format(3 / 2 * np.sqrt(k)))

# check measures of convergence
dEdchi_max = system.check_density_convergence('dEdchi')
mu_minus_dEdn_max = system.check_density_convergence('euler')
print('\nConvergence check:')
print('Max |𝛿E/𝛿χ| = {:.4g}'.format(dEdchi_max))
print('Max |µ-𝛿E/𝛿n| = {:.4g}'.format(mu_minus_dEdn_max))

This results in

Starting density optimization
  Iter      E [eV]      dE [eV]       Max |𝛿E/𝛿χ|       Max |µ-𝛿E/𝛿n|
   0     13613.510241      0            22.3543            999.713
   1     6558.513125     -7055          14.9388            60922.2
   2     3187.559139    -3370.95        10.6157            9693.09
   3     1626.061603    -1561.5         20.9591            10707.9
   4      889.606308    -736.455        6.05023            36347.1
   5      548.774621    -340.832        4.41879            16886.9
   6      376.605141    -172.169        3.18798            30169.6
   7      286.048990    -90.5562        13.3591          8.54989e+06
   8      230.137857    -55.9111        5.08921            14889.6
   9      193.892397    -36.2455        4.23115            21310.5
   10     169.331505    -24.5609        3.55797            10083.9
   11     152.787931    -16.5436        1.00027             269847
   12     143.654847    -9.13308        2.72997            56765.4
   13     135.496019    -8.15883        1.54196            91377.3
   14     132.008561    -3.48746         1.3482          2.23426e+06
   15     130.441073    -1.56749        0.514763            407990
   16     129.689363    -0.75171        1.57811             129319
   17     129.298829   -0.390534        0.993079           55204.9
   18     129.166347   -0.132481        0.241797            36499
   19     129.113693   -0.0526547       0.114015            932106
   20     129.094428   -0.0192641       0.104126            462777
   21     129.083642   -0.0107861       0.192647            134088
   22     129.078364  -0.00527784      0.0367585            204466
   23     129.076561  -0.00180324      0.0659604            221888
   24     129.075609  -0.000952299     0.0215513           67498.4
   25     129.075250  -0.000358949     0.0398467            559964
   26     129.075041  -0.000208968     0.0111308            471584
   27     129.074965  -7.56737e-05     0.00604273           307790
   28     129.074947  -1.8423e-05      0.00204187           104320
   29     129.074942  -5.12278e-06     0.00102419           459483
   30     129.074940  -1.90488e-06    0.000324086        1.26339e+06
   31     129.074939  -7.6968e-07     0.000166146        1.66163e+06
   32     129.074939  -3.65745e-07    0.000154676        1.66135e+07
   33     129.074939  -1.66777e-07    9.27869e-05           223466
   34     129.074939  -2.56686e-08    0.000148262           340946
   35     129.074939  -2.22304e-08     0.00010925        1.36216e+06
   36     129.074939  -2.2209e-08      0.00016803           592491
Density optimization successfully converged in 36 step(s)

Optimized energy = 4.74341650 Ha
Expected energy = 4.74341649 Ha

Convergence check:
Max |𝛿E/𝛿χ| = 0.0007487
Max |µ-𝛿E/𝛿n| = 5.925e+05

Note how the Max \(|\mu - \delta E / \delta n|\) measure of convergence becomes very large at convergence eventhough the optimized energy agrees with the theoretical one and Max \(|\delta E / \delta \chi|\) measure of convergence is small. This is a result of the vacuum regions with low densities leading to divergences in the von Weizsaecker term.

On this page
Show Source