Skip to content

Commit 4aa3954

Browse files
bastonerobastonero
committed
Add quasi-harmonic approximation scripts and notebooks
1 parent 6eb1ce3 commit 4aa3954

4 files changed

Lines changed: 888 additions & 0 deletions

File tree

‎Examples/sscha_and_aiida/model.ipynb‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,44 @@
3838
"})"
3939
]
4040
},
41+
{
42+
"cell_type": "code",
43+
"execution_count": null,
44+
"metadata": {},
45+
"outputs": [],
46+
"source": [
47+
"import numpy as np\n",
48+
"\n",
49+
"from aiida import load_profile\n",
50+
"from aiida.orm import *\n",
51+
"\n",
52+
"from qe_tools import CONSTANTS as C\n",
53+
"\n",
54+
"from ase.io import write, read\n",
55+
"from ase import units\n",
56+
"from ase.calculators.singlepoint import SinglePointCalculator\n",
57+
"\n",
58+
"from flare.atoms import FLARE_Atoms\n",
59+
"from flare.learners.utils import is_std_in_bound, get_env_indices\n",
60+
"from flare.bffs.sgp.calculator import SGP_Calculator\n",
61+
"from flare.bffs.sgp._C_flare import Structure\n",
62+
"\n",
63+
"import matplotlib.pyplot as plt\n",
64+
"\n",
65+
"load_profile()\n",
66+
"\n",
67+
"\n",
68+
"plt.rcParams.update({\n",
69+
" 'text.usetex': False,\n",
70+
" # 'text.latex.preamble': r'\\usepackage{sansmath} \\sansmath',\n",
71+
" 'pdf.fonttype':42,\n",
72+
" 'font.family':'sans-serif',\n",
73+
" 'font.sans-serif':'Arial',\n",
74+
" 'font.size':14,\n",
75+
" 'mathtext.fontset': 'stixsans',\n",
76+
"})"
77+
]
78+
},
4179
{
4280
"cell_type": "code",
4381
"execution_count": 2,

‎Examples/sscha_and_aiida/qha.ipynb‎

Lines changed: 535 additions & 0 deletions
Large diffs are not rendered by default.

‎Examples/sscha_and_aiida/qha.py‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
"""Script to run QHA with FLARE calculator."""
2+
import numpy as np
3+
4+
from ase import units, Atoms
5+
from ase.io import read
6+
from ase.optimize import BFGS
7+
from ase.constraints import ExpCellFilter
8+
9+
from phonopy import Phonopy
10+
from phonopy.structure.atoms import PhonopyAtoms
11+
12+
from flare.bffs.sgp.calculator import SGP_Calculator
13+
from flare.atoms import FLARE_Atoms
14+
15+
structure_filepath = 'Si.pwi'
16+
model_filepath = 'model.json'
17+
initial_vcrelax = True
18+
fmax = 0.001
19+
scale_i = 0.8
20+
scale_f = 1.2
21+
scale_step = 12
22+
supercell_matrix = [4, 4, 4]
23+
distance = 0.01
24+
t_min = 0
25+
t_max = 300
26+
symmetrize = True
27+
primitive_matrix = None
28+
mesh = 100
29+
unitcell0 = read(structure_filepath)
30+
31+
flare_calc, _ = SGP_Calculator.from_file(model_filepath)
32+
unitcell0.calc = flare_calc
33+
34+
def run() -> None:
35+
"""Run QHA calculation."""
36+
# Initial optimization
37+
print("Running initial geometry optimization...", flush=True)
38+
if initial_vcrelax:
39+
print("Initial volume: ", unitcell0.get_volume(), flush=True)
40+
unitcell = geometry_optimization(
41+
unitcell0,
42+
vcrelax=initial_vcrelax,
43+
fmax=fmax,
44+
)
45+
if initial_vcrelax:
46+
print("Final volume: ", unitcell.get_volume(), flush=True)
47+
48+
# Run equation of state
49+
print("Running equation of state...", flush=True)
50+
all_scaled_atoms, energies = run_eos(unitcell)
51+
52+
print("Running phonons for each volume...")
53+
start_index = -energies.index(min(energies))
54+
for scaled_atoms in all_scaled_atoms:
55+
phonons(scaled_atoms, filename=f'./thermal_properties.yaml-{start_index}')
56+
start_index += 1
57+
print("Run completed")
58+
59+
def run_eos(unitcell: Atoms) -> tuple[tuple[Atoms], list]:
60+
"""Run equation of states."""
61+
energies, volumes = [], []
62+
all_scaled_atoms = []
63+
64+
scale_factors = np.linspace(scale_i, scale_f, scale_step)
65+
66+
for scale_factor in scale_factors:
67+
scaled_atoms = scale_and_relax(unitcell, scale_factor, fmax=fmax)
68+
all_scaled_atoms.append(scaled_atoms)
69+
energies.append(scaled_atoms.get_potential_energy())
70+
volumes.append(scaled_atoms.get_volume())
71+
72+
np.savetxt('./e-v.dat', np.array([volumes, energies]).T)
73+
74+
return all_scaled_atoms, energies
75+
76+
def phonons(atoms: Atoms, filename: str = None) -> None:
77+
"""Run phonons using Phonopy."""
78+
unitcell = PhonopyAtoms(
79+
symbols=atoms.get_chemical_symbols(),
80+
numbers=atoms.get_atomic_numbers(),
81+
scaled_positions=atoms.get_scaled_positions(),
82+
cell=atoms.get_cell(),
83+
)
84+
85+
ph = Phonopy(
86+
unitcell=unitcell,
87+
primitive_matrix=primitive_matrix,
88+
supercell_matrix=supercell_matrix,
89+
)
90+
91+
ph.generate_displacements(distance=distance)
92+
supercells = ph.supercells_with_displacements
93+
94+
sets_of_forces = []
95+
for supercell in supercells:
96+
cell, scaled_positions, numbers = supercell.totuple()
97+
supercell_atoms = Atoms(
98+
cell=cell,
99+
scaled_positions=scaled_positions,
100+
numbers=numbers,
101+
calculator=flare_calc,
102+
)
103+
sets_of_forces.append(supercell_atoms.get_forces())
104+
105+
ph.forces = sets_of_forces
106+
ph.produce_force_constants()
107+
108+
if symmetrize:
109+
ph.symmetrize_force_constants()
110+
ph.symmetrize_force_constants_by_space_group()
111+
112+
ph.run_mesh(mesh=mesh)
113+
ph.run_thermal_properties(t_min=t_min, t_max=t_max)
114+
ph.thermal_properties.write_yaml(filename)
115+
116+
117+
def geometry_optimization(atoms: Atoms, vcrelax: bool = False, fmax: float = 0.001) -> Atoms:
118+
"""Optimize geometry of the given atoms."""
119+
if vcrelax:
120+
ecf = ExpCellFilter(atoms, scalar_pressure=0)
121+
optimizer = BFGS(ecf)
122+
else:
123+
optimizer = BFGS(atoms=atoms)
124+
125+
optimizer.run(fmax=fmax)
126+
127+
if vcrelax:
128+
return optimizer.atoms.atoms
129+
return optimizer.atoms
130+
131+
132+
def scale_and_relax(atoms: Atoms, scale_factor: float, fmax: float = 0.001) -> Atoms:
133+
"""Scale and relax the atoms of an ASE Atoms object."""
134+
ase = atoms.copy()
135+
ase.calc = atoms.calc
136+
ase.set_cell(ase.get_cell() * float(scale_factor) ** (1 / 3), scale_atoms=True)
137+
138+
return geometry_optimization(ase, vcrelax=False, fmax=fmax)
139+
140+
141+
if __name__ == "__main__":
142+
run()
143+

‎Examples/sscha_and_aiida/qha_flare.py‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
"""Script to run QHA with FLARE calculator."""
2+
import numpy as np
3+
4+
from ase import units, Atoms
5+
from ase.io import read
6+
from ase.optimize import BFGS
7+
from ase.constraints import ExpCellFilter
8+
9+
from phonopy import Phonopy
10+
from phonopy.structure.atoms import PhonopyAtoms
11+
12+
from flare.bffs.sgp.calculator import SGP_Calculator
13+
from flare.atoms import FLARE_Atoms
14+
15+
16+
class QHA:
17+
"""Run QHA with FLARE calculator."""
18+
19+
def __init__(
20+
self,
21+
structure_filepath = 'Si.pwi',
22+
model_filepath = 'model.json',
23+
initial_vcrelax = False,
24+
fmax = 0.001,
25+
unitcell = None,
26+
scale_i = 0.94,
27+
scale_f = 1.06,
28+
scale_step = 12,
29+
supercell_matrix = [2, 2, 2],
30+
distance = 0.01,
31+
t_min = 0,
32+
t_max = 300,
33+
symmetrize = False,
34+
primitive_matrix = None,
35+
mesh = 100,
36+
):
37+
"""Constructor of the class."""
38+
self.initial_vcrelax = initial_vcrelax
39+
self.fmax = fmax
40+
self.unitcell = FLARE_Atoms.from_ase_atoms(read(structure_filepath))
41+
self.scale_i = scale_i
42+
self.scale_f = scale_f
43+
self.scale_step = scale_step
44+
self.supercell_matrix = supercell_matrix
45+
self.distance = distance
46+
self.t_min = t_min
47+
self.t_max = t_max
48+
self.symmetrize = symmetrize
49+
self.primitive_matrix = primitive_matrix
50+
self.mesh = mesh
51+
52+
flare_calc, _ = SGP_Calculator.from_file(model_filepath)
53+
self.unitcell.calc = flare_calc
54+
55+
def run(self) -> None:
56+
"""Run QHA calculation."""
57+
print(self.unitcell.get_potential_energy(), flush=True)
58+
# Initial optimization
59+
print("Running initial geometry optimization...", flush=True)
60+
if self.initial_vcrelax:
61+
print("Initial volume: ", self.unitcell.get_volume(), flush=True)
62+
self.unitcell = geometry_optimization(
63+
self.unitcell,
64+
vcrelax=self.initial_vcrelax,
65+
fmax=self.fmax,
66+
)
67+
if self.initial_vcrelax:
68+
print("Final volume: ", self.unitcell.get_volume(), flush=True)
69+
70+
# Run equation of state
71+
print("Running equation of state...", flush=True)
72+
all_scaled_atoms, energies = self.run_eos()
73+
74+
print("Running phonons for each volume...")
75+
start_index = -energies.index(min(energies))
76+
for scaled_atoms in all_scaled_atoms:
77+
self.phonons(scaled_atoms, filename=f'./thermal_properties.yaml-{start_index}')
78+
start_index += 1
79+
80+
print("Run completed")
81+
82+
def run_eos(self) -> tuple[tuple[Atoms], list]:
83+
"""Run equation of states."""
84+
energies, volumes = [], []
85+
all_scaled_atoms = []
86+
87+
scale_factors = np.linspace(self.scale_i, self.scale_f, self.scale_step)
88+
89+
for scale_factor in scale_factors:
90+
scaled_atoms = scale_and_relax(self.unitcell, scale_factor, fmax=self.fmax)
91+
all_scaled_atoms.append(scaled_atoms)
92+
energies.append(scaled_atoms.get_potential_energy())
93+
volumes.append(scaled_atoms.get_volume())
94+
95+
np.savetxt('./e-v.dat', np.array([volumes, energies]).T)
96+
97+
return all_scaled_atoms, energies
98+
99+
def phonons(self, atoms: Atoms, filename: str = None) -> None:
100+
"""Run phonons using Phonopy."""
101+
unitcell = PhonopyAtoms(
102+
symbols=atoms.get_chemical_symbols(),
103+
numbers=atoms.get_atomic_numbers(),
104+
scaled_positions=atoms.get_scaled_positions(),
105+
cell=atoms.get_cell(),
106+
)
107+
108+
ph = Phonopy(
109+
unitcell=unitcell,
110+
primitive_matrix=self.primitive_matrix,
111+
supercell_matrix=self.supercell_matrix,
112+
)
113+
114+
ph.generate_displacements(distance=self.distance)
115+
supercells = ph.supercells_with_displacements
116+
117+
sets_of_forces = []
118+
for supercell in supercells:
119+
cell, scaled_positions, numbers = supercell.totuple()
120+
supercell_atoms = Atoms(
121+
cell=cell,
122+
scaled_positions=scaled_positions,
123+
numbers=numbers,
124+
calculator=self.flare_calc,
125+
)
126+
sets_of_forces.append(supercell_atoms.get_forces())
127+
128+
ph.forces = sets_of_forces
129+
ph.produce_force_constants()
130+
131+
if self.symmetrize:
132+
ph.symmetrize_force_constants()
133+
ph.symmetrize_force_constants_by_space_group()
134+
135+
ph.run_mesh(mesh=self.mesh)
136+
ph.run_thermal_properties(t_min=self.t_min, t_max=self.t_max)
137+
ph.thermal_properties.write_yaml(filename)
138+
139+
140+
def geometry_optimization(
141+
atoms: Atoms,
142+
vcrelax: bool = False,
143+
fmax: float = 0.001,
144+
) -> Atoms:
145+
"""Optimize geometry of the given atoms."""
146+
print("I am in opt", flush=True)
147+
print(atoms.calc)
148+
if vcrelax:
149+
ecf = ExpCellFilter(atoms, scalar_pressure=0)
150+
optimizer = BFGS(ecf)
151+
else:
152+
optimizer = BFGS(atoms=atoms)
153+
print("Running opt", flush=True)
154+
optimizer.run(fmax=fmax)
155+
print("I almost finished in opt", flush=True)
156+
157+
if vcrelax:
158+
return optimizer.atoms.atoms
159+
return optimizer.atoms
160+
161+
162+
def scale_and_relax(atoms: Atoms, scale_factor: float, fmax: float = 0.001) -> Atoms:
163+
"""Scale and relax the atoms of an ASE Atoms object."""
164+
ase = atoms.copy()
165+
ase.calc = atoms.calc
166+
ase.set_cell(ase.get_cell() * float(scale_factor) ** (1 / 3), scale_atoms=True)
167+
168+
return geometry_optimization(ase, vcrelax=False, fmax=fmax)
169+
170+
171+
QHA().run()
172+

0 commit comments

Comments
 (0)