Geometry Optimization#
Conventional Geometry Optimization#
The following code is an example of how a geometry optimization can be set up for body-centred cubic lithium (bcc-Li). This example is meant to illustrate how one can perform different types of geometry optimizations, including ones where
the lattice is fixed but the ions can move,
the lattice can deform but the fractional ionic coordinates are fixed, and
both the lattice and ions can be changed
import torch
from system import System
from functionals import IonIon, IonElectron, Hartree, WangTeterStyleFunctional, PerdewBurkeErnzerhof
# create system and compute ground state energy
box_len = 3.48
box_vecs = box_len * torch.eye(3, dtype=torch.double)
shape = System.ecut2shape(800, box_vecs)
ions = [['Li', 'li.gga.recpot', box_len * torch.tensor([[0, 0, 0], [0.5, 0.5, 0.5]]).double()]]
WTexp = WangTeterStyleFunctional((5 / 6, 5 / 6, lambda x: torch.exp(x)))
terms = [IonIon, IonElectron, Hartree, WTexp.forward, PerdewBurkeErnzerhof]
system = System(box_vecs, shape, ions, terms, units='a')
system.optimize_density(1e-10)
energy = system.energy('eV')
print('Initial Energy = {:.4f} eV per atom'.format(energy / system.ion_count()))
# OPTIMIZE IONIC POSITIONS / MINIMIZE FORCES
# peturb ions
print('Perturbing ions ...')
system.place_ions(box_len * torch.tensor([[0.0, 0.1, 0.0], [0.6, 0.4, 0.6]], dtype=torch.double), units='a')
system.optimize_density(1e-10)
print('Perturbed energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
# restore optimal ionic positions by minimizing forces, keeping lattice fixed
print('Performing force minimization ...')
system.optimize_geometry(stol=None, ftol=1e-4, g_method='LBFGSlinesearch', g_verbose=True)
print('Optimized Energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
# OPTIMIZE LATTICE / MINIMIZE STRESS
# predict relaxed energy by fitting to murnaghan equation
print('\nPerforming EOS fit for equilibrium volume ...')
params, err = system.eos_fit(N=5)
relaxed_energy = system.ion_count() * params[2]
print('Equilibrium energy = {:.4f} eV per atom'.format(relaxed_energy / system.ion_count()))
# distort lattice
print('Deforming lattice ...')
tm = torch.tensor([[0.94, -0.03, 0.05],
[-0.03, 0.98, 0.04],
[0.05, 0.04, 1.05]], dtype=torch.double)
system.set_lattice(torch.matmul(tm, system.lattice_vectors('a').T).T, units='a')
system.optimize_density(1e-10)
print('Perturbed energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
# relax by minimizing stress, keeping box coordinates of ions fixed
print('Performing stress minimization ...')
system.optimize_geometry(ftol=None, stol=1e-4, g_method='LBFGSlinesearch', g_verbose=True)
print('Optimized Energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
# OPTIMIZE GEOMETRY / MINIMIZE FORCES AND STRESS
print('\nPerturbing overall geometry ...')
# peturb ions and distort lattice
tm = torch.tensor([[0.94, -0.03, 0.05],
[-0.03, 0.98, 0.04],
[0.05, 0.04, 1.05]], dtype=torch.double)
system.place_ions(torch.matmul(tm, system.cartesian_ionic_coordinates('a').T).T, units='a')
system.set_lattice(torch.matmul(tm, system.lattice_vectors('a').T).T, units='a')
system.optimize_density(1e-10)
print('Perturbed energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
# restore optimal geometry by minimizing forces and stress
print('Performing geometry optimization ...')
system.optimize_geometry(stol=1e-4, ftol=1e-4, g_method='LBFGSlinesearch', g_verbose=True)
print('Optimized Energy = {:.4f} eV per atom'.format(system.energy('eV') / system.ion_count()))
This code results in the following output.
Initial Energy = -7.3578 eV per atom
Perturbing ions ...
Perturbed energy = -7.1217 eV per atom
Performing force minimization ...
Iter E [eV per atom] dE [eV per atom] Max Force [eV/Å] Max Stress [eV/ų]
0 -7.121696 0 0.842513 0.0185083
1 -7.310906 -0.189211 0.386264 0.00338438
2 -7.345332 -0.0344256 0.200588 0.00365126
3 -7.354303 -0.00897128 0.106939 0.00374702
4 -7.356813 -0.00251014 0.0572883 0.00379373
5 -7.357531 -0.000717481 0.0307186 0.00380226
6 -7.357738 -0.000206621 0.0164658 0.00380478
7 -7.357797 -5.91924e-05 0.00882966 0.00380545
8 -7.357814 -1.7049e-05 0.00473394 0.00380565
9 -7.357819 -4.90418e-06 0.00256241 0.00380571
10 -7.357820 -1.42075e-06 0.0013795 0.00380573
11 -7.357821 -4.08516e-07 0.000735595 0.00380574
12 -7.357821 -1.16475e-07 0.000396698 0.00380574
13 -7.357821 -3.37541e-08 0.000208924 0.00380575
14 -7.357821 -9.54252e-09 0.000110661 0.00380575
15 -7.357821 -4.41454e-09 1.77384e-05 0.00380575
16 -7.357821 -4.93685e-11 1.77835e-05 0.00380575
17 -7.357821 -3.02247e-11 1.7895e-05 0.00380575
Geometry optimization successfully converged in 17 step(s)
Optimized Energy = -7.3578 eV per atom
Performing EOS fit for equilibrium volume ...
Equilibrium energy = -7.3595 eV per atom
Deforming lattice ...
Perturbed energy = -7.3351 eV per atom
Performing stress minimization ...
Iter E [eV per atom] dE [eV per atom] Max Force [eV/Å] Max Stress [eV/ų]
0 -7.335086 0 0.00012742 0.00944689
1 -7.358675 -0.0235888 8.01787e-05 0.0015873
2 -7.358963 -0.000288282 5.3544e-05 0.00115537
3 -7.359343 -0.000380017 5.01887e-05 0.000458943
4 -7.359427 -8.35168e-05 4.60166e-05 0.000255862
5 -7.359440 -1.35338e-05 4.27197e-05 0.000363761
6 -7.359472 -3.224e-05 4.11179e-05 0.000176818
7 -7.359483 -1.10674e-05 3.99874e-05 8.06105e-05
8 -7.359487 -3.86039e-06 3.92089e-05 4.99504e-05
9 -7.359489 -2.03904e-06 3.85918e-05 3.05821e-05
Geometry optimization successfully converged in 9 step(s)
Optimized Energy = -7.3595 eV per atom
Perturbing overall geometry ...
Perturbed energy = -7.3068 eV per atom
Performing geometry optimization ...
Iter E [eV per atom] dE [eV per atom] Max Force [eV/Å] Max Stress [eV/ų]
0 -7.306823 0 0.396871 0.0116692
1 -7.335214 -0.0283912 0.26266 0.00706001
2 -7.353728 -0.0185138 0.126866 0.00332661
3 -7.357243 -0.00351547 0.0623917 0.00269824
4 -7.358423 -0.00118024 0.0434388 0.000865097
5 -7.358744 -0.000321098 0.0213561 0.000882875
6 -7.359002 -0.000257791 0.0111404 0.000413472
7 -7.359038 -3.54953e-05 0.00549463 0.000416077
8 -7.359083 -4.53939e-05 0.00307135 0.000412545
9 -7.359114 -3.08373e-05 0.00150737 0.000655856
10 -7.359122 -7.70214e-06 0.000644913 0.000701684
11 -7.359159 -3.71719e-05 0.0012657 0.000468603
12 -7.359179 -2.04144e-05 0.000537675 0.000531482
13 -7.359198 -1.87877e-05 0.000983285 0.000464925
14 -7.359246 -4.82697e-05 0.000741491 0.000537872
15 -7.359276 -2.935e-05 0.000590132 0.000664748
16 -7.359330 -5.46165e-05 0.000638868 0.000601527
17 -7.359360 -2.92485e-05 0.00941012 0.000724566
18 -7.359400 -4.03731e-05 0.00473564 0.000584351
19 -7.359438 -3.81269e-05 0.00231378 0.000216816
20 -7.359450 -1.2374e-05 0.000602998 0.00028038
21 -7.359474 -2.3748e-05 0.000323331 0.00011232
22 -7.359479 -4.69937e-06 0.000132017 0.000117203
23 -7.359484 -4.73172e-06 0.000899799 0.000158323
24 -7.359486 -2.43998e-06 0.000292098 0.00015559
25 -7.359486 -1.38117e-07 0.000260749 8.92468e-05
26 -7.359489 -2.7876e-06 0.000147261 6.82452e-05
27 -7.359491 -1.57216e-06 0.000107321 4.29506e-05
28 -7.359488 2.63341e-06 8.10592e-05 4.74667e-05
29 -7.359490 -2.3244e-06 5.92544e-05 3.53178e-05
30 -7.359492 -2.08458e-06 4.42381e-05 1.96232e-05
Geometry optimization successfully converged in 30 step(s)
Optimized Energy = -7.3595 eV per atom
Parameterized Geometry Optimization#
The following code is an example of how a parameterized geometry optimization can be set up for hexagonal close-packed magnesium (hcp-Mg).
import numpy as np
import torch
from system import System
from functionals import IonIon, IonElectron, Hartree, WangTeterStyleFunctional, PerdewBurkeErnzerhof
# use GPU if available else use CPU
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
params = torch.tensor([24 / System.A_per_b**3, 1.5], dtype=torch.double, device=device).requires_grad_()
print('Initial Guess: Volume per atom = {:.5f} ų, c/a = {:.5f}'
.format(params[0].item() * System.A_per_b**3, params[1].item()))
# define the lattice vectors and fractional ionic coordinates as a function of the parameters
def parameterized_geometry(params):
vol_per_atom, c_over_a = params
a = ((2 * torch.abs(vol_per_atom)) / (np.sqrt(3) / 2 * c_over_a)).pow(1 / 3)
box_vecs = torch.tensor([[1.0, 0.0, 0.0],
[-0.5, np.sqrt(3) / 2, 0.0],
[0.0, 0.0, 0.0]], dtype=torch.double, device=device)
box_vecs[2, 2] = torch.abs(c_over_a)
box_vecs = a * box_vecs
frac_ion_coords = torch.tensor([[1 / 3, 2 / 3, 3 / 4],
[2 / 3, 1 / 3, 1 / 4]], dtype=torch.double, device=device)
return box_vecs, frac_ion_coords
box_vecs, frac_ion_coords = parameterized_geometry(params)
WTexp = WangTeterStyleFunctional((5 / 6, 5 / 6, lambda x: torch.exp(x)))
# required for GPU usage with functionals that inherit from the KineticFunctional class
WTexp.set_device(device)
terms = [IonIon, IonElectron, Hartree, WTexp.forward, PerdewBurkeErnzerhof]
# construct the system object
ions = [['Mg', 'mg.gga.recpot', frac_ion_coords.detach()]]
# lattice vectors must be in angstroms for ecut2shape
shape = System.ecut2shape(2000, box_vecs.detach() * System.A_per_b)
system = System(box_vecs, shape, ions, terms, units='b', coord_type='fractional', device=device)
# define a print statement to track how the parameters evolve over the optimization
def param_string(params):
return '{:.5f} {:.5f}'.format(params[0].item() * System.A_per_b**3, params[1].item())
system.optimize_parameterized_geometry(params, parameterized_geometry, g_method='LBFGSlinesearch',
g_verbose=True, param_string=param_string, ftol=1e-3, stol=1e-3)
print('Optimized Results: Volume per atom = {:.5f} ų, c/a = {:.5f}\n'
.format(params[0].item() * System.A_per_b**3, params[1].item()))
This code results in the following output.
Initial Guess: Volume per atom = 24.00000 ų, c/a = 1.50000
Iter E [eV per atom] dE [eV per atom] Max Force [eV/Å] Max Stress [eV/ų] Params
0 -24.198576 0 3.87304e-06 0.0138319 24.00000 1.50000
1 -24.208906 -0.0103299 4.73783e-06 0.01135 23.99641 1.55999
2 -24.214733 -0.00582784 4.90708e-06 0.00547138 23.48824 1.59311
3 -24.216084 -0.00135011 5.02009e-06 0.00247967 23.28204 1.61115
4 -24.216418 -0.000333999 5.09427e-06 0.00116472 23.21050 1.62080
5 -24.216488 -7.04247e-05 5.12923e-06 0.000593536 23.16811 1.62585
6 -24.216517 -2.87005e-05 5.15291e-06 0.000389594 23.15198 1.62861
7 -24.216529 -1.2568e-05 5.17015e-06 0.000188783 23.15116 1.63011
Geometry optimization successfully converged in 7 step(s)
Optimized Results: Volume per atom = 23.15116 ų, c/a = 1.63011