Workflows#
To demonstrate the workflows implemented in the atomistics package, the LAMMPS molecular
dynamics simulation code is used in the following demonstrations. Still the same workflows can also be used with other
simulation codes:
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = evaluate_with_lammpslib(
task_dict={},
potential_dataframe=potential_dataframe,
)
/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/calculators/__init__.py:63: UserWarning: calc_static_with_qe(), evaluate_with_qe() and optimize_positions_and_volume_with_qe() are not available as the import of the module named 'pwtools' failed.
raise_warning(module_list=quantum_espresso_function, import_error=e)
The interatomic potential for Aluminium from Mishin named 1999--Mishin-Y--Al--LAMMPS--ipr1 is used in the evaluation
with LAMMPS evaluate_with_lammpslib().
Elastic Matrix#
The elastic constants and elastic moduli can be calculated using the ElasticMatrixWorkflow:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import (
analyse_results_for_elastic_matrix,
get_tasks_for_elastic_matrix,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
task_dict, sym_dict = get_tasks_for_elastic_matrix(
structure=bulk("Al", cubic=True),
num_of_point=5,
eps_range=0.005,
sqrt_eta=True,
)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
fit_dict, sym_dict = analyse_results_for_elastic_matrix(output_dict=result_dict, sym_dict=sym_dict)
print(fit_dict)
{'elastic_matrix': array([[114.10393022, 60.51098897, 60.51098897, 0. ,
0. , 0. ],
[ 60.51098897, 114.10393022, 60.51098897, 0. ,
0. , 0. ],
[ 60.51098897, 60.51098897, 114.10393022, 0. ,
0. , 0. ],
[ 0. , 0. , 0. , 31.67539618,
0. , 0. ],
[ 0. , 0. , 0. , 0. ,
31.67539618, 0. ],
[ 0. , 0. , 0. , 0. ,
0. , 31.67539618]]), 'elastic_matrix_inverse': array([[ 0.01385713, -0.00480204, -0.00480204, 0. , 0. ,
0. ],
[-0.00480204, 0.01385713, -0.00480204, 0. , 0. ,
0. ],
[-0.00480204, -0.00480204, 0.01385713, 0. , 0. ,
0. ],
[ 0. , 0. , 0. , 0.03157024, 0. ,
0. ],
[ 0. , 0. , 0. , 0. , 0.03157024,
0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0.03157024]]), 'bulkmodul_voigt': 78.37530271895345, 'bulkmodul_reuss': 78.37530271895342, 'bulkmodul_hill': 78.37530271895344, 'shearmodul_voigt': 29.72382596007553, 'shearmodul_reuss': 29.525101479037495, 'shearmodul_hill': 29.624463719556513, 'youngsmodul_voigt': 79.16385035531343, 'youngsmodul_reuss': 78.69363099993072, 'youngsmodul_hill': 78.92882891555043, 'poissonsratio_voigt': 0.3316564708332767, 'poissonsratio_reuss': 0.3326564018044502, 'poissonsratio_hill': 0.3321562486791239, 'AVR': 0.3354060396152373, 'elastic_matrix_eigval': EigResult(eigenvalues=array([ 53.59294125, 235.12590816, 53.59294125, 31.67539618,
31.67539618, 31.67539618]), eigenvectors=array([[-0.81649658, 0.57735027, 0.33985586, 0. , 0. ,
0. ],
[ 0.40824829, 0.57735027, -0.81286921, 0. , 0. ,
0. ],
[ 0.40824829, 0.57735027, 0.47301336, 0. , 0. ,
0. ],
[ 0. , 0. , 0. , 1. , 0. ,
0. ],
[ 0. , 0. , 0. , 0. , 1. ,
0. ],
[ 0. , 0. , 0. , 0. , 0. ,
1. ]]))}
The ElasticMatrixWorkflow takes an ase.atoms.Atoms object as structure input as well as the number of points
num_of_point for each compression direction. Depending on the symmetry of the input structure the number of
calculations required to calculate the elastic matrix changes. The compression and elongation range is defined by the
eps_range parameter. Furthermore, sqrt_eta and fit_order describe how the change in energy over compression and
elongation is fitted to calculate the resulting pressure.
Energy Volume Curve#
The EnergyVolumeCurveWorkflow can be used to calculate the equilibrium properties: equilibrium volume, equilibrium
energy, equilibrium bulk modulus and the pressure derivative of the equilibrium bulk modulus.
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import (
analyse_results_for_energy_volume_curve,
get_tasks_for_energy_volume_curve,
get_thermal_properties_for_energy_volume_curve,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
structure = bulk("Al", cubic=True)
task_dict = get_tasks_for_energy_volume_curve(
structure=structure,
num_points=11,
vol_range=0.05,
axes=("x", "y", "z"),
)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
fit_dict = analyse_results_for_energy_volume_curve(
output_dict=result_dict,
task_dict=task_dict,
fit_type="polynomial",
fit_order=3,
)
print(fit_dict)
{'b_prime_eq': 1.279502459079921, 'bulkmodul_eq': 77.72501358958125, 'volume_eq': 66.43019853103964, 'energy_eq': -13.43996804374383, 'fit_dict': {'fit_type': 'polynomial', 'least_square_error': 3.225313797039607e-10, 'poly_fit': array([-4.17645808e-05, 1.19746500e-02, -1.03803906e+00, 1.49168639e+01]), 'fit_order': 3}, 'energy': [-13.398169481534445, -13.413389552957456, -13.425112589013958, -13.433411420804067, -13.438357630783006, -13.439999952539933, -13.438383476946305, -13.433607982916406, -13.425774537190858, -13.414961805921427, -13.401233093668836], 'volume': [63.10861874999998, 63.77291999999998, 64.43722124999998, 65.1015225, 65.76582375000004, 66.43012500000002, 67.09442624999994, 67.75872750000002, 68.42302874999999, 69.08732999999997, 69.75163125000002]}
The input parameters for the EnergyVolumeCurveWorkflow in addition to the ase.atoms.Atoms object defined
as structure are:
num_pointsthe number of strains to calculate energies and volumes.fit_typethe type of the fit which should be used to calculate the equilibrium properties. This can either be apolynomialfit or a specific equation of state like the Birch equation (birch), the Birch-Murnaghan equation (birchmurnaghan) the Murnaghan equation (murnaghan), the Pourier Tarnatola eqaution (pouriertarantola) or the Vinet equation (vinet).fit_orderfor thepolynomialfit type the order of the polynomial can be set, for the other fit types this parameter is ignored.vol_rangespecifies the amount of compression and elongation to be applied relative to the absolute volume.axesspecifies the axes which are compressed, typically a uniform compression is applied.strainsspecifies the strains directly rather than deriving them from the range of volume compressionvol_range.
Beyond calculating the equilibrium properties the EnergyVolumeCurveWorkflow can also be used to calculate the thermal
properties using the Moruzzi, V. L. et al. model:
tp_dict = get_thermal_properties_for_energy_volume_curve(
fit_dict=fit_dict,
masses=structure.get_masses(),
t_min=1,
t_max=1500,
t_step=50,
temperatures=None,
constant_volume=False,
)
print(tp_dict)
{'temperatures': array([ 1, 51, 101, 151, 201, 251, 301, 351, 401, 451, 501,
551, 601, 651, 701, 751, 801, 851, 901, 951, 1001, 1051,
1101, 1151, 1201, 1251, 1301, 1351, 1401, 1451, 1501]), 'volumes': array([66.48459155, 66.48492729, 66.48841343, 66.49613572, 66.50654263,
66.51846055, 66.53126421, 66.5446199 , 66.55833931, 66.57230985,
66.58646057, 66.6007448 , 66.61513063, 66.6295956 , 66.64412341,
66.65870199, 66.6733222 , 66.68797701, 66.70266093, 66.71736958,
66.73209946, 66.74684773, 66.76161205, 66.77639048, 66.79118142,
66.8059835 , 66.82079558, 66.83561668, 66.85044595, 66.86528267,
66.88012622]), 'free_energy': array([ 0.18879418, 0.18840183, 0.18352524, 0.16909367, 0.1440755 ,
0.10931095, 0.06593656, 0.01498215, -0.04269081, -0.1063728 ,
-0.1754776 , -0.24951635, -0.328077 , -0.41080851, -0.49740877,
-0.58761537, -0.68119851, -0.77795536, -0.87770572, -0.98028844,
-1.08555864, -1.19338539, -1.3036498 , -1.41624343, -1.53106703,
-1.6480294 , -1.76704645, -1.88804043, -2.01093923, -2.13567578,
-2.26218757]), 'entropy': array([ 0.75685476, 5.08219061, 18.62461551, 38.05446425,
57.66932288, 75.37710503, 90.99476551, 104.78762775,
117.06473009, 128.09164492, 138.08127287, 147.20167192,
155.5857919 , 163.33970924, 170.54896549, 177.28330935,
183.6002256 , 189.54757242, 195.16556895, 200.48830823,
205.5449212 , 210.36048156, 214.95671658, 219.35257074,
223.56465686, 227.60762031, 231.49443545, 235.23664865,
238.84457906, 242.32748553, 244.04031817]), 'heat_capacity': array([8.65067171e-02, 9.11255799e+00, 3.33019964e+01, 5.89575081e+01,
7.50185080e+01, 8.36468610e+01, 8.85256734e+01, 9.15055757e+01,
9.34491088e+01, 9.47846079e+01, 9.57412353e+01, 9.64498999e+01,
9.69896043e+01, 9.74102601e+01, 9.77446368e+01, 9.80149634e+01,
9.82367471e+01, 9.84210719e+01, 9.85760297e+01, 9.87076399e+01,
9.88204550e+01, 9.89179695e+01, 9.90029019e+01, 9.90773926e+01,
9.91431454e+01, 9.92015302e+01, 9.92536586e+01, 9.93004401e+01,
9.93426247e+01, nan, nan])}
Or alternatively directly calculate the thermal expansion:
thermal_properties_dict = get_thermal_properties_for_energy_volume_curve(
fit_dict=fit_dict,
masses=structure.get_masses(),
t_min=1,
t_max=1500,
t_step=50,
constant_volume=False,
output_keys=["temperatures", "volumes"],
)
temperatures, volumes = (
thermal_properties_dict["temperatures"],
thermal_properties_dict["volumes"],
)
The Moruzzi, V. L. et al. model is a quantum mechanical approximation, so the equilibrium volume at 0K is not the same as the equilibrium volume calculated by fitting the equation of state.
Molecular Dynamics#
Just like the structure optimization also the molecular dynamics calculation can either be implemented inside the
simulation code or in the atomistics package. The latter has the advantage that it is the same implementation for all
different simulation codes, while the prior has the advantage that it is usually faster and computationally more efficient.
Implemented in Simulation Code#
The LAMMPS simulation code implements a wide range of different simulation workflows, this
includes molecular dynamics. In the atomistics package these can be directly accessed via the python interface.
Langevin Thermostat#
The Langevin thermostat is currently the only thermostat which is available as both a stand-alone python interface and an integrated interface inside the LAMMPS simulation code. The latter is introduced here:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_langevin_with_lammpslib,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = calc_molecular_dynamics_langevin_with_lammpslib(
structure=bulk("Al", cubic=True).repeat([10, 10, 10]),
potential_dataframe=potential_dataframe,
Tstart=100,
Tstop=100,
Tdamp=0.1,
run=100,
thermo=10,
timestep=0.001,
seed=4928459,
dist="gaussian",
output_keys=(
"positions",
"cell",
"forces",
"temperature",
"energy_pot",
"energy_tot",
"pressure",
"velocities",
),
)
In addition to the typical LAMMPS input parameters like the atomistic structure structure as ase.atoms.Atoms object
and the pandas.DataFrame for the interatomic potential potential_dataframe are:
Tstartstart temperatureTstopend temperatureTdamptemperature damping parameterrunnumber of molecular dynamics steps to be executed during one temperature stepthermorefresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular dynamics steps.timesteptime step - typically 1fs defined as0.001.seedrandom seed for the molecular dynamicsdistinitial velocity distributionlmpLammps library instance aspylammpsmpi.LammpsASELibraryobjectoutputthe output quantities which are extracted from the molecular dynamics simulation
Nose Hoover Thermostat#
Canonical ensemble (nvt) - volume and temperature constraints molecular dynamics:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_nvt_with_lammpslib,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = calc_molecular_dynamics_nvt_with_lammpslib(
structure=bulk("Al", cubic=True).repeat([10, 10, 10]),
potential_dataframe=potential_dataframe,
Tstart=100,
Tstop=100,
Tdamp=0.1,
run=100,
thermo=10,
timestep=0.001,
seed=4928459,
dist="gaussian",
output_keys=(
"positions",
"cell",
"forces",
"temperature",
"energy_pot",
"energy_tot",
"pressure",
),
)
In addition to the typical LAMMPS input parameters like the atomistic structure structure as ase.atoms.Atoms object
and the pandas.DataFrame for the interatomic potential potential_dataframe are:
Tstartstart temperatureTstopend temperatureTdamptemperature damping parameterrunnumber of molecular dynamics steps to be executed during one temperature stepthermorefresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular dynamics steps.timesteptime step - typically 1fs defined as0.001.seedrandom seed for the molecular dynamicsdistinitial velocity distributionlmpLammps library instance aspylammpsmpi.LammpsASELibraryobjectoutputthe output quantities which are extracted from the molecular dynamics simulation
Isothermal-isobaric ensemble (npt) - pressure and temperature constraints molecular dynamics:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_npt_with_lammpslib,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = calc_molecular_dynamics_npt_with_lammpslib(
structure=bulk("Al", cubic=True).repeat([10, 10, 10]),
potential_dataframe=potential_dataframe,
Tstart=100,
Tstop=100,
Tdamp=0.1,
run=100,
thermo=100,
timestep=0.001,
Pstart=0.0,
Pstop=0.0,
Pdamp=1.0,
seed=4928459,
dist="gaussian",
output_keys=(
"positions",
"cell",
"forces",
"temperature",
"energy_pot",
"energy_tot",
"pressure",
),
)
The input parameters for the isothermal-isobaric ensemble (npt) are the same as for the canonical ensemble (nvt) plus:
Pstartstart pressurePstopend pressurePdamppressure damping parameter
Isenthalpic ensemble (nph) - pressure and helmholtz-energy constraints molecular dynamics:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_nph_with_lammpslib,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = calc_molecular_dynamics_nph_with_lammpslib(
structure=bulk("Al", cubic=True).repeat([10, 10, 10]),
potential_dataframe=potential_dataframe,
run=100,
thermo=100,
timestep=0.001,
Tstart=100,
Pstart=0.0,
Pstop=0.0,
Pdamp=1.0,
seed=4928459,
dist="gaussian",
output_keys=(
"positions",
"cell",
"forces",
"temperature",
"energy_pot",
"energy_tot",
"pressure",
),
)
Thermal Expansion#
One example of a molecular dynamics calculation with the LAMMPS simulation code is the calculation of the thermal expansion:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_thermal_expansion_with_lammpslib,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
temperatures_md, volumes_md = calc_molecular_dynamics_thermal_expansion_with_lammpslib(
structure=bulk("Al", cubic=True).repeat([10, 10, 10]),
potential_dataframe=potential_dataframe,
Tstart=100,
Tstop=1000,
Tstep=100,
Tdamp=0.1,
run=100,
thermo=100,
timestep=0.001,
Pstart=0.0,
Pstop=0.0,
Pdamp=1.0,
seed=4928459,
dist="gaussian",
lmp=None,
)
100%|██████████| 10/10 [00:05<00:00, 1.79it/s]
In addition to the typical LAMMPS input parameters like the atomistic structure structure as ase.atoms.Atoms object
and the pandas.DataFrame for the interatomic potential potential_dataframe are:
Tstartstart temperatureTstopend temperatureTsteptemperature stepTdamptemperature damping parameterrunnumber of molecular dynamics steps to be executed during one temperature stepthermorefresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular dynamics steps.timesteptime step - typically 1fs defined as0.001.Pstartstart pressurePstopend pressurePdamppressure damping parameterseedrandom seed for the molecular dynamicsdistinitial velocity distributionlmpLammps library instance aspylammpsmpi.LammpsASELibraryobject
These input parameters are based on the LAMMPS fix nvt/npt, you can read more about the specific implementation on the
LAMMPS website.
Phonons from Molecular Dynamics#
The softening of the phonon modes is calculated for Silicon using the Tersoff interatomic potential which is available via the NIST potentials repository. Silicon is chosen based on its diamond crystal lattice which requires less calculation than the face centered cubic (fcc) crystal of Aluminium. The simulation workflow consists of three distinct steps:
Starting with the optimization of the equilibrium structure.
Followed by the calculation of the 0K phonon spectrum.
Finally, the finite temperature phonon spectrum is calculated using molecular dynamics.
The finite temperature phonon spectrum is calculated using the DynaPhoPy
package, which is integrated inside the atomistics package. As a prerequisite the dependencies, imported and the bulk
silicon diamond structure is created and the Tersoff interatomic potential is loaded:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_phonons_with_lammpslib,
evaluate_with_lammpslib,
)
from atomistics.workflows import (
optimize_positions_and_volume,
get_tasks_for_harmonic_approximation,
analyse_results_for_harmonic_approximation,
)
from dynaphopy import Quasiparticle
import pandas
import spglib
structure_bulk = bulk("Si", cubic=True)
potential_dataframe = get_potential_by_name(
potential_name="1988--Tersoff-J--Si-c--LAMMPS--ipr1", resource_path="static/lammps"
)
The first step is optimizing the Silicon diamond structure to match the lattice specifications implemented in the Tersoff interatomic potential:
task_dict = optimize_positions_and_volume(structure=structure_bulk)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
structure_ase = result_dict["structure_with_optimized_positions_and_volume"]
As a second step the 0K phonons are calculated using the PhonopyWorkflow which is explained in more detail below in
the section on Phonons.
cell = (
structure_ase.cell.array,
structure_ase.get_scaled_positions(),
structure_ase.numbers,
)
primitive_matrix = spglib.standardize_cell(cell=cell, to_primitive=True)[
0
] / structure_ase.get_volume() ** (1 / 3)
task_dict, phonopy_obj = get_tasks_for_harmonic_approximation(
structure=structure_ase,
interaction_range=10,
displacement=0.01,
primitive_matrix=primitive_matrix,
number_of_snapshots=None,
)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
phonopy_dict = analyse_results_for_harmonic_approximation(
phonopy=phonopy_obj,
output_dict=result_dict,
dos_mesh=20,
number_of_snapshots=None,
)
The calcualtion of the finite temperature phonons starts by computing the molecular dynamics trajectory using the
calc_molecular_dynamics_phonons_with_lammpslib() function. This function is internally linked to DynaPhoPy
to return an dynaphopy.dynamics.Dynamics object:
trajectory = calc_molecular_dynamics_phonons_with_lammpslib(
structure_ase=structure_ase,
potential_dataframe=potential_dataframe,
force_constants=phonopy_obj.force_constants,
phonopy_unitcell=phonopy_obj.unitcell,
phonopy_primitive_matrix=phonopy_obj.primitive_matrix,
phonopy_supercell_matrix=phonopy_obj.supercell_matrix,
total_time=2, # ps
time_step=0.001, # ps
relaxation_time=5, # ps
silent=True,
supercell=[2, 2, 2],
memmap=False,
velocity_only=True,
temperature=600,
)
When a total of 2 picoseconds is selected to compute the finite temperature phonons with a timestep of 1 femto second then this results in a total of 2000 molecular dynamics steps. While more molecular dynamics steps result in more precise predictions they also require more computational resources.
The postprocessing is executed using the DynaPhoPy package:
calculation = Quasiparticle(trajectory)
calculation.select_power_spectra_algorithm(2) # select FFT algorithm
calculation.get_renormalized_phonon_dispersion_bands()
renormalized_force_constants = (
calculation.get_renormalized_force_constants().get_array()
)
renormalized_force_constants
Using 2000 steps
Using Fast Fourier transform (Numpy) function
set frequency range: 0.0 - 21.200000000000003
Q-point: 1 / 32 [ 0.00000 0.00000 0.00000 ]
Harmonic frequencies (THz):
[-4.12629462e-06 -4.12313400e-06 -4.11976059e-06 1.60678991e+01
1.60678991e+01 1.60678991e+01]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
MD cell size relation: [2 2 2]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Peak # 1
----------------------------------------------
Width 0.471003 THz
Position 0.032516 THz
Area (<K>) (Lorentzian) 0.000000 eV
Area (<K>) (Total) 0.000000 eV
<|dQ/dt|^2> 0.000000 eV
Occupation number -0.500000
Fit temperature nan K
Base line -0.000000 eV * ps
Maximum height 0.000000 eV * ps
Fitting global error 814564757159.949097
Frequency shift 0.032520 THz
Peak # 2
----------------------------------------------
Width 0.471003 THz
Position 0.032516 THz
Area (<K>) (Lorentzian) 0.000000 eV
Area (<K>) (Total) 0.000000 eV
<|dQ/dt|^2> 0.000000 eV
Occupation number -0.500000
Fit temperature nan K
Base line -0.000000 eV * ps
Maximum height 0.000000 eV * ps
Fitting global error 814564757159.949097
Frequency shift 0.032520 THz
Peak # 3
----------------------------------------------
Width 0.471003 THz
Position 0.032516 THz
Area (<K>) (Lorentzian) 0.000000 eV
Area (<K>) (Total) 0.000000 eV
<|dQ/dt|^2> 0.000000 eV
Occupation number -0.500000
Fit temperature nan K
Base line -0.000000 eV * ps
Maximum height 0.000000 eV * ps
Fitting global error 814564757159.949097
Frequency shift 0.032520 THz
Peak # 4
----------------------------------------------
Width 0.790596 THz
Position 15.562379 THz
Area (<K>) (Lorentzian) 0.014512 eV
Area (<K>) (Total) 0.014213 eV
<|dQ/dt|^2> 0.029024 eV
Occupation number 2.333456
Fit temperature 333.284829 K
Base line 0.000007 eV * ps
Maximum height 0.011686 eV * ps
Fitting global error 0.032669
Frequency shift -0.505520 THz
Peak # 5
----------------------------------------------
Width 0.790596 THz
Position 15.562379 THz
Area (<K>) (Lorentzian) 0.014512 eV
Area (<K>) (Total) 0.014213 eV
<|dQ/dt|^2> 0.029024 eV
Occupation number 2.333456
Fit temperature 333.284829 K
Base line 0.000007 eV * ps
Maximum height 0.011686 eV * ps
Fitting global error 0.032669
Frequency shift -0.505520 THz
Peak # 6
----------------------------------------------
Width 0.790596 THz
Position 15.562379 THz
Area (<K>) (Lorentzian) 0.014512 eV
Area (<K>) (Total) 0.014213 eV
<|dQ/dt|^2> 0.029024 eV
Occupation number 2.333456
Fit temperature 333.284829 K
Base line 0.000007 eV * ps
Maximum height 0.011686 eV * ps
Fitting global error 0.032669
Frequency shift -0.505520 THz
Fixing gamma point 0 frequencies
Q-point: 2 / 32 [ 0.00000 0.25000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Peak # 1
----------------------------------------------
Width 0.521044 THz
Position 4.512483 THz
Area (<K>) (Lorentzian) 0.018118 eV
Area (<K>) (Total) 0.016810 eV
<|dQ/dt|^2> 0.036236 eV
Occupation number 11.699901
Fit temperature 420.263285 K
Base line -0.000042 eV * ps
Maximum height 0.022137 eV * ps
Fitting global error 0.016824
Frequency shift -0.151490 THz
Peak # 2
----------------------------------------------
Width 0.521044 THz
Position 4.512483 THz
Area (<K>) (Lorentzian) 0.018118 eV
Area (<K>) (Total) 0.016810 eV
<|dQ/dt|^2> 0.036236 eV
Occupation number 11.699901
Fit temperature 420.263285 K
Base line -0.000042 eV * ps
Maximum height 0.022137 eV * ps
Fitting global error 0.016824
Frequency shift -0.151490 THz
Peak # 3
----------------------------------------------
Width 0.884643 THz
Position 6.802090 THz
Area (<K>) (Lorentzian) 0.034381 eV
Area (<K>) (Total) 0.042075 eV
<|dQ/dt|^2> 0.068762 eV
Occupation number 14.858240
Fit temperature 797.669962 K
Base line 0.000413 eV * ps
Maximum height 0.024742 eV * ps
Fitting global error 0.034947
Frequency shift -0.096079 THz
Peak # 4
----------------------------------------------
Width 0.816224 THz
Position 14.710909 THz
Area (<K>) (Lorentzian) 0.050561 eV
Area (<K>) (Total) 0.056117 eV
<|dQ/dt|^2> 0.101122 eV
Occupation number 9.943370
Fit temperature 1172.575803 K
Base line 0.000331 eV * ps
Maximum height 0.039435 eV * ps
Fitting global error 0.029537
Frequency shift -0.459579 THz
Peak # 5
----------------------------------------------
Width 0.848633 THz
Position 15.068935 THz
Area (<K>) (Lorentzian) 0.021045 eV
Area (<K>) (Total) 0.019693 eV
<|dQ/dt|^2> 0.042091 eV
Occupation number 3.743663
Fit temperature 486.177110 K
Base line -0.000033 eV * ps
Maximum height 0.015788 eV * ps
Fitting global error 0.028720
Frequency shift -0.486744 THz
Peak # 6
----------------------------------------------
Width 0.848633 THz
Position 15.068935 THz
Area (<K>) (Lorentzian) 0.021045 eV
Area (<K>) (Total) 0.019693 eV
<|dQ/dt|^2> 0.042091 eV
Occupation number 3.743663
Fit temperature 486.177110 K
Base line -0.000033 eV * ps
Maximum height 0.015788 eV * ps
Fitting global error 0.028720
Frequency shift -0.486744 THz
Q-point: 3 / 32 [ 0.00000 0.50000 0.50000 ]
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Peak # 1
----------------------------------------------
Width 0.581916 THz
Position 6.560644 THz
Area (<K>) (Lorentzian) 0.025435 eV
Area (<K>) (Total) 0.027461 eV
<|dQ/dt|^2> 0.050870 eV
Occupation number 11.280117
Fit temperature 589.967377 K
Base line 0.000120 eV * ps
Maximum height 0.027826 eV * ps
Fitting global error 0.025135
Frequency shift -0.334692 THz
Peak # 2
----------------------------------------------
Width 0.581916 THz
Position 6.560644 THz
Area (<K>) (Lorentzian) 0.025435 eV
Area (<K>) (Total) 0.027461 eV
<|dQ/dt|^2> 0.050870 eV
Occupation number 11.280117
Fit temperature 589.967377 K
Base line 0.000120 eV * ps
Maximum height 0.027826 eV * ps
Fitting global error 0.025135
Frequency shift -0.334692 THz
Peak # 3
----------------------------------------------
Width 0.603434 THz
Position 11.933607 THz
Area (<K>) (Lorentzian) 0.030449 eV
Area (<K>) (Total) 0.034566 eV
<|dQ/dt|^2> 0.060899 eV
Occupation number 7.253040
Fit temperature 705.721741 K
Base line 0.000221 eV * ps
Maximum height 0.032124 eV * ps
Fitting global error 0.023757
Frequency shift -0.258183 THz
Peak # 4
----------------------------------------------
Width 0.603434 THz
Position 11.933607 THz
Area (<K>) (Lorentzian) 0.030449 eV
Area (<K>) (Total) 0.034566 eV
<|dQ/dt|^2> 0.060899 eV
Occupation number 7.253040
Fit temperature 705.721741 K
Base line 0.000221 eV * ps
Maximum height 0.032124 eV * ps
Fitting global error 0.023757
Frequency shift -0.258183 THz
Peak # 5
----------------------------------------------
Width 0.633146 THz
Position 14.446483 THz
Area (<K>) (Lorentzian) 0.042340 eV
Area (<K>) (Total) 0.047530 eV
<|dQ/dt|^2> 0.084680 eV
Occupation number 8.405362
Fit temperature 981.634126 K
Base line 0.000288 eV * ps
Maximum height 0.042572 eV * ps
Fitting global error 0.015345
Frequency shift -0.444472 THz
Peak # 6
----------------------------------------------
Width 0.633146 THz
Position 14.446483 THz
Area (<K>) (Lorentzian) 0.042340 eV
Area (<K>) (Total) 0.047530 eV
<|dQ/dt|^2> 0.084680 eV
Occupation number 8.405362
Fit temperature 981.634126 K
Base line 0.000288 eV * ps
Maximum height 0.042572 eV * ps
Fitting global error 0.015345
Frequency shift -0.444472 THz
Q-point: 4 / 32 [ 0.00000 0.75000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Peak # 1
----------------------------------------------
Width 0.514465 THz
Position 4.501657 THz
Area (<K>) (Lorentzian) 0.024087 eV
Area (<K>) (Total) 0.021361 eV
<|dQ/dt|^2> 0.048173 eV
Occupation number 15.758007
Fit temperature 558.850479 K
Base line -0.000102 eV * ps
Maximum height 0.029806 eV * ps
Fitting global error 0.014263
Frequency shift -0.162316 THz
Peak # 2
----------------------------------------------
Width 0.514465 THz
Position 4.501657 THz
Area (<K>) (Lorentzian) 0.024087 eV
Area (<K>) (Total) 0.021361 eV
<|dQ/dt|^2> 0.048173 eV
Occupation number 15.758007
Fit temperature 558.850479 K
Base line -0.000102 eV * ps
Maximum height 0.029806 eV * ps
Fitting global error 0.014263
Frequency shift -0.162316 THz
Peak # 3
----------------------------------------------
Width 0.840450 THz
Position 6.833587 THz
Area (<K>) (Lorentzian) 0.047750 eV
Area (<K>) (Total) 0.056965 eV
<|dQ/dt|^2> 0.095499 eV
Occupation number 20.731750
Fit temperature 1108.018836 K
Base line 0.000500 eV * ps
Maximum height 0.036169 eV * ps
Fitting global error 0.027488
Frequency shift -0.064582 THz
Peak # 4
----------------------------------------------
Width 0.864892 THz
Position 14.761576 THz
Area (<K>) (Lorentzian) 0.036911 eV
Area (<K>) (Total) 0.042399 eV
<|dQ/dt|^2> 0.073823 eV
Occupation number 7.097870
Fit temperature 855.439739 K
Base line 0.000312 eV * ps
Maximum height 0.027169 eV * ps
Fitting global error 0.033172
Frequency shift -0.408912 THz
Peak # 5
----------------------------------------------
Width 0.688214 THz
Position 15.046611 THz
Area (<K>) (Lorentzian) 0.029921 eV
Area (<K>) (Total) 0.029439 eV
<|dQ/dt|^2> 0.059841 eV
Occupation number 5.542238
Fit temperature 692.843575 K
Base line 0.000013 eV * ps
Maximum height 0.027678 eV * ps
Fitting global error 0.016812
Frequency shift -0.509068 THz
Peak # 6
----------------------------------------------
Width 0.688214 THz
Position 15.046611 THz
Area (<K>) (Lorentzian) 0.029921 eV
Area (<K>) (Total) 0.029439 eV
<|dQ/dt|^2> 0.059841 eV
Occupation number 5.542238
Fit temperature 692.843575 K
Base line 0.000013 eV * ps
Maximum height 0.027678 eV * ps
Fitting global error 0.016812
Frequency shift -0.509068 THz
Q-point: 5 / 32 [ 0.25000 0.00000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Skipped, equivalent to [0. 0.25 0.25]
Q-point: 6 / 32 [ 0.25000 0.25000 0.50000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Peak # 1
----------------------------------------------
Width 0.530326 THz
Position 4.473769 THz
Area (<K>) (Lorentzian) 0.045467 eV
Area (<K>) (Total) 0.044412 eV
<|dQ/dt|^2> 0.090934 eV
Occupation number 30.380653
Fit temperature 1055.152004 K
Base line 0.000002 eV * ps
Maximum height 0.054580 eV * ps
Fitting global error 0.012315
Frequency shift -0.194110 THz
Peak # 2
----------------------------------------------
Width 0.887612 THz
Position 6.692894 THz
Area (<K>) (Lorentzian) 0.080224 eV
Area (<K>) (Total) 0.092591 eV
<|dQ/dt|^2> 0.160448 eV
Occupation number 35.921329
Fit temperature 1861.810062 K
Base line 0.000700 eV * ps
Maximum height 0.057539 eV * ps
Fitting global error 0.020878
Frequency shift -0.268196 THz
Peak # 3
----------------------------------------------
Width 0.880545 THz
Position 8.797643 THz
Area (<K>) (Lorentzian) 0.042403 eV
Area (<K>) (Total) 0.050339 eV
<|dQ/dt|^2> 0.084807 eV
Occupation number 14.145261
Fit temperature 983.756386 K
Base line 0.000429 eV * ps
Maximum height 0.030657 eV * ps
Fitting global error 0.028868
Frequency shift -0.208203 THz
Peak # 4
----------------------------------------------
Width 0.759458 THz
Position 13.374413 THz
Area (<K>) (Lorentzian) 0.019163 eV
Area (<K>) (Total) 0.022883 eV
<|dQ/dt|^2> 0.038327 eV
Occupation number 3.853723
Fit temperature 442.800477 K
Base line 0.000197 eV * ps
Maximum height 0.016064 eV * ps
Fitting global error 0.039426
Frequency shift -0.350502 THz
Peak # 5
----------------------------------------------
Width 0.649011 THz
Position 14.982356 THz
Area (<K>) (Lorentzian) 0.024299 eV
Area (<K>) (Total) 0.024730 eV
<|dQ/dt|^2> 0.048598 eV
Occupation number 4.427993
Fit temperature 562.012770 K
Base line 0.000048 eV * ps
Maximum height 0.023835 eV * ps
Fitting global error 0.012789
Frequency shift -0.444090 THz
Peak # 6
----------------------------------------------
Width 0.788135 THz
Position 15.142930 THz
Area (<K>) (Lorentzian) 0.061707 eV
Area (<K>) (Total) 0.071316 eV
<|dQ/dt|^2> 0.123414 eV
Occupation number 11.881954
Fit temperature 1431.382786 K
Base line 0.000537 eV * ps
Maximum height 0.049844 eV * ps
Fitting global error 0.022823
Frequency shift -0.439826 THz
Q-point: 7 / 32 [ 0.25000 0.50000 0.75000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Peak # 1
----------------------------------------------
Width 0.906301 THz
Position 7.246475 THz
Area (<K>) (Lorentzian) 0.079723 eV
Area (<K>) (Total) 0.098689 eV
<|dQ/dt|^2> 0.159446 eV
Occupation number 32.928885
Fit temperature 1850.159723 K
Base line 0.001009 eV * ps
Maximum height 0.056001 eV * ps
Fitting global error 0.021638
Frequency shift -0.296021 THz
Peak # 2
----------------------------------------------
Width 0.906301 THz
Position 7.246475 THz
Area (<K>) (Lorentzian) 0.079723 eV
Area (<K>) (Total) 0.098689 eV
<|dQ/dt|^2> 0.159446 eV
Occupation number 32.928885
Fit temperature 1850.159723 K
Base line 0.001009 eV * ps
Maximum height 0.056001 eV * ps
Fitting global error 0.021638
Frequency shift -0.296021 THz
Peak # 3
----------------------------------------------
Width 0.640129 THz
Position 11.066449 THz
Area (<K>) (Lorentzian) 0.020915 eV
Area (<K>) (Total) 0.024452 eV
<|dQ/dt|^2> 0.041830 eV
Occupation number 5.242680
Fit temperature 484.188840 K
Base line 0.000186 eV * ps
Maximum height 0.020800 eV * ps
Fitting global error 0.025943
Frequency shift -0.283872 THz
Peak # 4
----------------------------------------------
Width 0.640129 THz
Position 11.066449 THz
Area (<K>) (Lorentzian) 0.020915 eV
Area (<K>) (Total) 0.024452 eV
<|dQ/dt|^2> 0.041830 eV
Occupation number 5.242680
Fit temperature 484.188840 K
Base line 0.000186 eV * ps
Maximum height 0.020800 eV * ps
Fitting global error 0.025943
Frequency shift -0.283872 THz
Peak # 5
----------------------------------------------
Width 0.819240 THz
Position 14.812285 THz
Area (<K>) (Lorentzian) 0.040303 eV
Area (<K>) (Total) 0.044734 eV
<|dQ/dt|^2> 0.080605 eV
Occupation number 7.767513
Fit temperature 934.242317 K
Base line 0.000265 eV * ps
Maximum height 0.031319 eV * ps
Fitting global error 0.031010
Frequency shift -0.426053 THz
Peak # 6
----------------------------------------------
Width 0.819240 THz
Position 14.812285 THz
Area (<K>) (Lorentzian) 0.040303 eV
Area (<K>) (Total) 0.044734 eV
<|dQ/dt|^2> 0.080605 eV
Occupation number 7.767513
Fit temperature 934.242317 K
Base line 0.000265 eV * ps
Maximum height 0.031319 eV * ps
Fitting global error 0.031010
Frequency shift -0.426053 THz
Q-point: 8 / 32 [ 0.25000 0.75000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Peak # 1
----------------------------------------------
Width 0.525422 THz
Position 4.481425 THz
Area (<K>) (Lorentzian) 0.035386 eV
Area (<K>) (Total) 0.033182 eV
<|dQ/dt|^2> 0.070771 eV
Occupation number 23.492554
Fit temperature 821.150247 K
Base line -0.000064 eV * ps
Maximum height 0.042875 eV * ps
Fitting global error 0.013007
Frequency shift -0.186454 THz
Peak # 2
----------------------------------------------
Width 0.897054 THz
Position 6.708112 THz
Area (<K>) (Lorentzian) 0.055174 eV
Area (<K>) (Total) 0.065058 eV
<|dQ/dt|^2> 0.110348 eV
Occupation number 24.491888
Fit temperature 1280.366541 K
Base line 0.000548 eV * ps
Maximum height 0.039156 eV * ps
Fitting global error 0.026906
Frequency shift -0.252979 THz
Peak # 3
----------------------------------------------
Width 0.830756 THz
Position 8.809692 THz
Area (<K>) (Lorentzian) 0.011457 eV
Area (<K>) (Total) 0.014099 eV
<|dQ/dt|^2> 0.022914 eV
Occupation number 3.451648
Fit temperature 264.483182 K
Base line 0.000139 eV * ps
Maximum height 0.008780 eV * ps
Fitting global error 0.062783
Frequency shift -0.196154 THz
Peak # 4
----------------------------------------------
Width 0.844979 THz
Position 13.308378 THz
Area (<K>) (Lorentzian) 0.031581 eV
Area (<K>) (Total) 0.038344 eV
<|dQ/dt|^2> 0.063163 eV
Occupation number 6.710551
Fit temperature 731.794210 K
Base line 0.000359 eV * ps
Maximum height 0.023794 eV * ps
Fitting global error 0.034462
Frequency shift -0.416538 THz
Peak # 5
----------------------------------------------
Width 0.667099 THz
Position 14.951007 THz
Area (<K>) (Lorentzian) 0.033338 eV
Area (<K>) (Total) 0.034961 eV
<|dQ/dt|^2> 0.066676 eV
Occupation number 6.275402
Fit temperature 772.339381 K
Base line 0.000114 eV * ps
Maximum height 0.031815 eV * ps
Fitting global error 0.015770
Frequency shift -0.475439 THz
Peak # 6
----------------------------------------------
Width 0.945956 THz
Position 15.121793 THz
Area (<K>) (Lorentzian) 0.019154 eV
Area (<K>) (Total) 0.019823 eV
<|dQ/dt|^2> 0.038308 eV
Occupation number 3.348779
Fit temperature 442.036286 K
Base line 0.000063 eV * ps
Maximum height 0.012891 eV * ps
Fitting global error 0.031504
Frequency shift -0.460962 THz
Q-point: 9 / 32 [ 0.50000 0.00000 0.50000 ]
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Skipped, equivalent to [0. 0.5 0.5]
Q-point: 10 / 32 [ 0.50000 0.25000 0.75000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Skipped, equivalent to [0.25 0.5 0.75]
Q-point: 11 / 32 [ 0.50000 0.50000 0.00000 ]
Harmonic frequencies (THz):
[ 6.89533567 6.89533567 12.19179039 12.19179039 14.89095524 14.89095524]
Skipped, equivalent to [0. 0.5 0.5]
Q-point: 12 / 32 [ 0.50000 0.75000 0.25000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Skipped, equivalent to [0.25 0.5 0.75]
Q-point: 13 / 32 [ 0.75000 0.00000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Skipped, equivalent to [0. 0.75 0.75]
Q-point: 14 / 32 [ 0.75000 0.25000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0.75 0. ]
Q-point: 15 / 32 [ 0.75000 0.50000 0.25000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Skipped, equivalent to [0.25 0.5 0.75]
Q-point: 16 / 32 [ 0.75000 0.75000 0.50000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0.75 0. ]
Q-point: 17 / 32 [ 0.25000 0.25000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Skipped, equivalent to [0.25 0. 0.25]
Q-point: 18 / 32 [ 0.25000 0.50000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0.25 0.5 ]
Q-point: 19 / 32 [ 0.25000 0.75000 0.50000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Skipped, equivalent to [0.25 0.5 0.75]
Q-point: 20 / 32 [ 0.25000 0.00000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.75 0.75 0.5 ]
Q-point: 21 / 32 [ 0.50000 0.25000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0.5 0.25]
Q-point: 22 / 32 [ 0.50000 0.50000 0.50000 ]
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Calculating phonon projection power spectra
Projecting into phonon mode
Projecting into wave vector
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Projecting into phonon mode
Projecting into wave vector
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Power spectrum resolution requested unavailable, using maximum: 0.500000 THz
If you need higher resolution increase the number of data
FFT: [##############################] 100.00% Done...
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Peak # 1
----------------------------------------------
Width 0.530025 THz
Position 4.474583 THz
Area (<K>) (Lorentzian) 0.041318 eV
Area (<K>) (Total) 0.038746 eV
<|dQ/dt|^2> 0.082636 eV
Occupation number 27.557693
Fit temperature 958.851771 K
Base line -0.000074 eV * ps
Maximum height 0.049628 eV * ps
Fitting global error 0.013256
Frequency shift -0.193296 THz
Peak # 2
----------------------------------------------
Width 0.530025 THz
Position 4.474583 THz
Area (<K>) (Lorentzian) 0.041318 eV
Area (<K>) (Total) 0.038746 eV
<|dQ/dt|^2> 0.082636 eV
Occupation number 27.557693
Fit temperature 958.851771 K
Base line -0.000074 eV * ps
Maximum height 0.049628 eV * ps
Fitting global error 0.013256
Frequency shift -0.193296 THz
Peak # 3
----------------------------------------------
Width 0.541104 THz
Position 10.987927 THz
Area (<K>) (Lorentzian) 0.053102 eV
Area (<K>) (Total) 0.050474 eV
<|dQ/dt|^2> 0.106205 eV
Occupation number 14.184593
Fit temperature 1231.977168 K
Base line -0.000083 eV * ps
Maximum height 0.062476 eV * ps
Fitting global error 0.009049
Frequency shift -0.323287 THz
Peak # 4
----------------------------------------------
Width 0.822189 THz
Position 12.813795 THz
Area (<K>) (Lorentzian) 0.015671 eV
Area (<K>) (Total) 0.018460 eV
<|dQ/dt|^2> 0.031342 eV
Occupation number 3.216043
Fit temperature 361.501217 K
Base line 0.000151 eV * ps
Maximum height 0.012134 eV * ps
Fitting global error 0.051115
Frequency shift -0.341043 THz
Peak # 5
----------------------------------------------
Width 0.667520 THz
Position 14.943763 THz
Area (<K>) (Lorentzian) 0.054206 eV
Area (<K>) (Total) 0.059370 eV
<|dQ/dt|^2> 0.108411 eV
Occupation number 10.521681
Fit temperature 1257.194693 K
Base line 0.000305 eV * ps
Maximum height 0.051696 eV * ps
Fitting global error 0.013639
Frequency shift -0.482683 THz
Peak # 6
----------------------------------------------
Width 0.667520 THz
Position 14.943763 THz
Area (<K>) (Lorentzian) 0.054206 eV
Area (<K>) (Total) 0.059370 eV
<|dQ/dt|^2> 0.108411 eV
Occupation number 10.521681
Fit temperature 1257.194693 K
Base line 0.000305 eV * ps
Maximum height 0.051696 eV * ps
Fitting global error 0.013639
Frequency shift -0.482683 THz
Q-point: 23 / 32 [ 0.50000 0.75000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0. 0.75]
Q-point: 24 / 32 [ 0.50000 0.00000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Skipped, equivalent to [0.5 0.5 0.5]
Q-point: 25 / 32 [ 0.75000 0.25000 0.50000 ]
Harmonic frequencies (THz):
[ 7.54249562 7.54249562 11.3503204 11.3503204 15.23833788 15.23833788]
Skipped, equivalent to [0.25 0.5 0.75]
Q-point: 26 / 32 [ 0.75000 0.50000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0. 0.75]
Q-point: 27 / 32 [ 0.75000 0.75000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66397327 4.66397327 6.89816884 15.17048811 15.55567884 15.55567884]
Skipped, equivalent to [0. 0.75 0.75]
Q-point: 28 / 32 [ 0.75000 0.00000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0. 0.75]
Q-point: 29 / 32 [ 0.00000 0.25000 0.75000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0. 0.75]
Q-point: 30 / 32 [ 0.00000 0.50000 0.00000 ]
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Skipped, equivalent to [0.5 0.5 0.5]
Q-point: 31 / 32 [ 0.00000 0.75000 0.25000 ]
Harmonic frequencies (THz):
[ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543]
Skipped, equivalent to [0.25 0. 0.75]
Q-point: 32 / 32 [ 0.00000 0.00000 0.50000 ]
Harmonic frequencies (THz):
[ 4.66787904 4.66787904 11.31121369 13.15483786 15.42644585 15.42644585]
Skipped, equivalent to [0.5 0.5 0.5]
array([[[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]],
[[-6.12113397e-03, -5.07266507e-04, -5.07266507e-04],
[-5.07266507e-04, 9.52223449e-03, 5.07266507e-04],
[-5.07266507e-04, 5.07266507e-04, 9.52223449e-03]],
[[ 9.52223449e-03, -5.07266507e-04, 5.07266507e-04],
[-5.07266507e-04, -6.12113397e-03, -5.07266507e-04],
[ 5.07266507e-04, -5.07266507e-04, 9.52223449e-03]],
...,
[[-5.80821903e-03, 1.30918037e-03, -2.34637025e-03],
[ 1.30918037e-03, 4.51532631e-03, 1.21434882e-02],
[-2.34637025e-03, 1.21434882e-02, 6.17293624e-03]],
[[ 4.51532631e-03, 1.30918037e-03, 1.21434882e-02],
[ 1.30918037e-03, -5.80821903e-03, -2.34637025e-03],
[ 1.21434882e-02, -2.34637025e-03, 6.17293624e-03]],
[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]]],
[[[-6.12113397e-03, -5.07266507e-04, -5.07266507e-04],
[-5.07266507e-04, 9.52223449e-03, 5.07266507e-04],
[-5.07266507e-04, 5.07266507e-04, 9.52223449e-03]],
[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]],
[[-6.96355984e-03, 5.07266507e-04, -5.07266507e-04],
[ 5.07266507e-04, -6.96355984e-03, -5.07266507e-04],
[-5.07266507e-04, -5.07266507e-04, 3.17249493e-03]],
...,
[[-3.48253365e+00, -2.36508256e+00, 2.36611975e+00],
[-2.36508256e+00, -3.48253365e+00, 2.36611975e+00],
[ 2.36611975e+00, 2.36611975e+00, -3.48419126e+00]],
[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]],
[[ 4.51532631e-03, 1.30918037e-03, 1.21434882e-02],
[ 1.30918037e-03, -5.80821903e-03, -2.34637025e-03],
[ 1.21434882e-02, -2.34637025e-03, 6.17293624e-03]]],
[[[ 9.52223449e-03, -5.07266507e-04, 5.07266507e-04],
[-5.07266507e-04, -6.12113397e-03, -5.07266507e-04],
[ 5.07266507e-04, -5.07266507e-04, 9.52223449e-03]],
[[-6.96355984e-03, 5.07266507e-04, -5.07266507e-04],
[ 5.07266507e-04, -6.96355984e-03, -5.07266507e-04],
[-5.07266507e-04, -5.07266507e-04, 3.17249493e-03]],
[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]],
...,
[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]],
[[-3.48253365e+00, -2.36508256e+00, 2.36611975e+00],
[-2.36508256e+00, -3.48253365e+00, 2.36611975e+00],
[ 2.36611975e+00, 2.36611975e+00, -3.48419126e+00]],
[[-5.80821903e-03, 1.30918037e-03, -2.34637025e-03],
[ 1.30918037e-03, 4.51532631e-03, 1.21434882e-02],
[-2.34637025e-03, 1.21434882e-02, 6.17293624e-03]]],
...,
[[[-5.80821903e-03, 1.30918037e-03, -2.34637025e-03],
[ 1.30918037e-03, 4.51532631e-03, 1.21434882e-02],
[-2.34637025e-03, 1.21434882e-02, 6.17293624e-03]],
[[-3.48253365e+00, -2.36508256e+00, 2.36611975e+00],
[-2.36508256e+00, -3.48253365e+00, 2.36611975e+00],
[ 2.36611975e+00, 2.36611975e+00, -3.48419126e+00]],
[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]],
...,
[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]],
[[-6.96355984e-03, 5.07266507e-04, -5.07266507e-04],
[ 5.07266507e-04, -6.96355984e-03, -5.07266507e-04],
[-5.07266507e-04, -5.07266507e-04, 3.17249493e-03]],
[[ 9.52223449e-03, -5.07266507e-04, 5.07266507e-04],
[-5.07266507e-04, -6.12113397e-03, -5.07266507e-04],
[ 5.07266507e-04, -5.07266507e-04, 9.52223449e-03]]],
[[[ 4.51532631e-03, 1.30918037e-03, 1.21434882e-02],
[ 1.30918037e-03, -5.80821903e-03, -2.34637025e-03],
[ 1.21434882e-02, -2.34637025e-03, 6.17293624e-03]],
[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]],
[[-3.48253365e+00, -2.36508256e+00, 2.36611975e+00],
[-2.36508256e+00, -3.48253365e+00, 2.36611975e+00],
[ 2.36611975e+00, 2.36611975e+00, -3.48419126e+00]],
...,
[[-6.96355984e-03, 5.07266507e-04, -5.07266507e-04],
[ 5.07266507e-04, -6.96355984e-03, -5.07266507e-04],
[-5.07266507e-04, -5.07266507e-04, 3.17249493e-03]],
[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]],
[[-6.12113397e-03, -5.07266507e-04, -5.07266507e-04],
[-5.07266507e-04, 9.52223449e-03, 5.07266507e-04],
[-5.07266507e-04, 5.07266507e-04, 9.52223449e-03]]],
[[[ 6.30333769e-03, -3.30497443e-05, 4.87343513e-03],
[-3.30497443e-05, 6.30333769e-03, 4.87343513e-03],
[ 4.87343513e-03, 4.87343513e-03, -1.18770702e-03]],
[[ 4.51532631e-03, 1.30918037e-03, 1.21434882e-02],
[ 1.30918037e-03, -5.80821903e-03, -2.34637025e-03],
[ 1.21434882e-02, -2.34637025e-03, 6.17293624e-03]],
[[-5.80821903e-03, 1.30918037e-03, -2.34637025e-03],
[ 1.30918037e-03, 4.51532631e-03, 1.21434882e-02],
[-2.34637025e-03, 1.21434882e-02, 6.17293624e-03]],
...,
[[ 9.52223449e-03, -5.07266507e-04, 5.07266507e-04],
[-5.07266507e-04, -6.12113397e-03, -5.07266507e-04],
[ 5.07266507e-04, -5.07266507e-04, 9.52223449e-03]],
[[-6.12113397e-03, -5.07266507e-04, -5.07266507e-04],
[-5.07266507e-04, 9.52223449e-03, 5.07266507e-04],
[-5.07266507e-04, 5.07266507e-04, 9.52223449e-03]],
[[ 1.48011729e+01, 5.07266507e-04, 5.07266507e-04],
[ 5.07266507e-04, 1.48011729e+01, 5.07266507e-04],
[ 5.07266507e-04, 5.07266507e-04, 1.48011729e+01]]]])
It calculates the re-normalized force constants which can then be used to calculate the finite temperature properties.
In addition the DynaPhoPy package can be used to directly compare the finite temperature phonon spectrum with the 0K phonon spectrum calulated with the finite displacement method:
calculation.plot_renormalized_phonon_dispersion_bands()
Langevin Thermostat#
In addition to the molecular dynamics implemented in the LAMMPS simulation code, the atomistics package also provides
the LangevinWorkflow which implements molecular dynamics independent of the specific simulation code.
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib_library_interface, get_potential_by_name
from atomistics.workflows import LangevinWorkflow
from pylammpsmpi import LammpsASELibrary
steps = 300
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
workflow = LangevinWorkflow(
structure=bulk("Al", cubic=True).repeat([2, 2, 2]),
temperature=1000.0,
overheat_fraction=2.0,
damping_timescale=100.0,
time_step=1,
)
lmp = LammpsASELibrary(
working_directory=None,
cores=1,
comm=None,
logger=None,
log_file=None,
library=None,
disable_log_file=True,
)
eng_pot_lst, eng_kin_lst = [], []
for i in range(steps):
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammpslib_library_interface(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
lmp=lmp,
)
eng_pot, eng_kin = workflow.analyse_structures(output_dict=result_dict)
eng_pot_lst.append(eng_pot)
eng_kin_lst.append(eng_kin)
lmp.close()
The advantage of this implementation is that the user can directly interact with the simulation between the individual
molecular dynamics simulation steps. This provides a lot of flexibility to prototype new simulation methods. The input
parameters of the LangevinWorkflow are:
structurethease.atoms.Atomsobject which is used as initial structure for the molecular dynamics calculationtemperaturethe temperature of the molecular dynamics calculation given in Kelvinoverheat_fractionthe over heating fraction of the Langevin thermostatdamping_timescalethe damping timescale of the Langevin thermostattime_stepthe time steps of the Langevin thermostat
Harmonic Approximation#
The harmonic approximation is implemented in two variations, once with constant volume and once including the volume expansion at finite temperature also known as quasi-harmonic approximation. Both of these are based on the phonopy package.
Phonons#
To calculate the phonons at a fixed volume the PhonopyWorkflow is used:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import (
get_band_structure,
get_hesse_matrix,
get_dynamical_matrix,
get_thermal_properties_for_harmonic_approximation,
get_tasks_for_harmonic_approximation,
analyse_results_for_harmonic_approximation,
plot_band_structure,
plot_dos,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
task_dict, phonopy_obj = get_tasks_for_harmonic_approximation(
structure=bulk("Al", cubic=True),
interaction_range=10,
displacement=0.01,
primitive_matrix=None,
number_of_snapshots=None,
)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
phonopy_dict = analyse_results_for_harmonic_approximation(
phonopy=phonopy_obj,
output_dict=result_dict,
dos_mesh=20,
number_of_snapshots=None,
)
The PhonopyWorkflow takes the following inputs:
structurethease.atoms.Atomsobject to calculate the phonon spectruminteraction_rangethe cutoff radius to consider for identifying the interaction between the atomsdisplacementdisplacement to calculate the forcesdos_meshmesh for the density of statesprimitive_matrixprimitive matrixnumber_of_snapshotsnumber of snapshots to calculate
In addition to the phonon properties, the PhonopyWorkflow also enables the calculation of thermal properties:
tp_dict = get_thermal_properties_for_harmonic_approximation(
phonopy=phonopy_obj,
t_min=1,
t_max=1500,
t_step=50,
temperatures=None,
cutoff_frequency=None,
pretend_real=False,
band_indices=None,
is_projection=False,
)
print(tp_dict)
{'temperatures': array([1.000e+00, 5.100e+01, 1.010e+02, 1.510e+02, 2.010e+02, 2.510e+02,
3.010e+02, 3.510e+02, 4.010e+02, 4.510e+02, 5.010e+02, 5.510e+02,
6.010e+02, 6.510e+02, 7.010e+02, 7.510e+02, 8.010e+02, 8.510e+02,
9.010e+02, 9.510e+02, 1.001e+03, 1.051e+03, 1.101e+03, 1.151e+03,
1.201e+03, 1.251e+03, 1.301e+03, 1.351e+03, 1.401e+03, 1.451e+03,
1.501e+03]), 'volumes': array([66.430125, 66.430125, 66.430125, 66.430125, 66.430125, 66.430125,
66.430125, 66.430125, 66.430125, 66.430125, 66.430125, 66.430125,
66.430125, 66.430125, 66.430125, 66.430125, 66.430125, 66.430125,
66.430125, 66.430125, 66.430125, 66.430125, 66.430125, 66.430125,
66.430125, 66.430125, 66.430125, 66.430125, 66.430125, 66.430125,
66.430125]), 'free_energy': array([ 0.14914132, 0.14837894, 0.13954171, 0.11738723, 0.08264779,
0.03712237, -0.01759836, -0.08025513, -0.14986079, -0.22563203,
-0.30693668, -0.39325592, -0.48415731, -0.57927552, -0.67829812,
-0.78095507, -0.88701079, -0.99625805, -1.10851315, -1.22361223,
-1.3414082 , -1.46176834, -1.58457228, -1.70971039, -1.8370824 ,
-1.96659625, -2.09816715, -2.23171671, -2.3671723 , -2.5044664 ,
-2.64353611]), 'entropy': array([1.10363871e-08, 5.98829810e+00, 2.96478195e+01, 5.54593816e+01,
7.80099308e+01, 9.71787932e+01, 1.13608521e+02, 1.27894607e+02,
1.40492150e+02, 1.51738264e+02, 1.61883985e+02, 1.71119149e+02,
1.79589851e+02, 1.87410480e+02, 1.94672040e+02, 2.01447985e+02,
2.07798389e+02, 2.13772961e+02, 2.19413270e+02, 2.24754417e+02,
2.29826293e+02, 2.34654555e+02, 2.39261386e+02, 2.43666089e+02,
2.47885561e+02, 2.51934678e+02, 2.55826598e+02, 2.59573021e+02,
2.63184393e+02, 2.66670075e+02, 2.70038493e+02]), 'heat_capacity': array([1.78544597e-07, 1.73410821e+01, 5.37349237e+01, 7.35976295e+01,
8.34733324e+01, 8.87978444e+01, 9.19287453e+01, 9.39060819e+01,
9.52277477e+01, 9.61520364e+01, 9.68225162e+01, 9.73237288e+01,
9.77079209e+01, 9.80087218e+01, 9.82485402e+01, 9.84427587e+01,
9.86022130e+01, 9.87347097e+01, 9.88459861e+01, 9.89403338e+01,
9.90210141e+01, 9.90905402e+01, 9.91508741e+01, 9.92035655e+01,
9.92498509e+01, 9.92907269e+01, 9.93270039e+01, 9.93593459e+01,
9.93883017e+01, 9.94143276e+01, 9.94378055e+01])}
The calculation of the thermal properties takes additional inputs:
t_minminimum temperaturet_maxmaximum temperaturet_steptemperature steptemperaturesalternative tot_min,t_maxandt_stepthe array of temperatures can be defined directlycutoff_frequencycutoff frequency to exclude the contributions of frequencies below a certain cut offpretend_realuse the absolute values of the phonon frequenciesband_indicesselect bands based on their indicesis_projectionmultiplies the squared eigenvectors - not recommended
Furthermore, also the dynamical matrix can be directly calculated with the PhonopyWorkflow:
mat = get_dynamical_matrix(phonopy=phonopy_obj, npoints=101)
mat
array([[ 1.72794621e-01, 6.42929783e-20, -6.42929783e-20,
-1.55150365e-18, 0.00000000e+00, -1.50515236e-19,
2.05737531e-18, -1.60732446e-20, 9.98736570e-02,
8.87128656e-18, -2.27745360e-34, -1.50493730e-19],
[ 6.42929783e-20, 1.40379905e-01, 1.92878935e-19,
-7.13795448e-36, 8.80589092e-18, -2.46493644e-33,
0.00000000e+00, 0.00000000e+00, -1.60732446e-20,
-1.11525914e-37, 8.80589092e-18, -4.46122155e-37],
[-6.42929783e-20, 1.92878935e-19, 1.72794621e-01,
-1.50493730e-19, 9.61616306e-34, 8.87128656e-18,
9.98736570e-02, 0.00000000e+00, -8.22950122e-18,
-1.50515236e-19, -1.47220311e-35, -1.55150365e-18],
[-1.55150365e-18, -7.13795448e-36, -1.50493730e-19,
1.72794621e-01, -1.92878935e-19, 6.42929783e-20,
8.87128656e-18, -2.29083727e-34, -1.50493730e-19,
8.22950122e-18, 0.00000000e+00, 9.98736570e-02],
[ 0.00000000e+00, 8.80589092e-18, 9.61616306e-34,
-1.92878935e-19, 1.40379905e-01, 6.02746672e-20,
2.51189080e-33, 8.80589092e-18, 1.78448862e-36,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[-1.50515236e-19, -2.46493644e-33, 8.87128656e-18,
6.42929783e-20, 6.02746672e-20, 1.72794621e-01,
-1.50515236e-19, 4.46122155e-37, -1.55150365e-18,
9.98736570e-02, 0.00000000e+00, -8.22950122e-18],
[ 2.05737531e-18, 0.00000000e+00, 9.98736570e-02,
8.87128656e-18, 2.51189080e-33, -1.50515236e-19,
1.72794621e-01, 1.92878935e-19, -6.42929783e-20,
-1.55150365e-18, 0.00000000e+00, -1.50493730e-19],
[-1.60732446e-20, 0.00000000e+00, 0.00000000e+00,
-2.29083727e-34, 8.80589092e-18, 4.46122155e-37,
1.92878935e-19, 1.40379905e-01, 6.02746672e-20,
-5.56657141e-48, 8.80589092e-18, 1.32252913e-33],
[ 9.98736570e-02, -1.60732446e-20, -8.22950122e-18,
-1.50493730e-19, 1.78448862e-36, -1.55150365e-18,
-6.42929783e-20, 6.02746672e-20, 1.72794621e-01,
-1.50515236e-19, -3.58804896e-33, 8.87128656e-18],
[ 8.87128656e-18, -1.11525914e-37, -1.50515236e-19,
8.22950122e-18, 0.00000000e+00, 9.98736570e-02,
-1.55150365e-18, -5.56657141e-48, -1.50515236e-19,
1.72794621e-01, 6.42929783e-20, -6.42929783e-20],
[-2.27745360e-34, 8.80589092e-18, -1.47220311e-35,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 8.80589092e-18, -3.58804896e-33,
6.42929783e-20, 1.40379905e-01, 1.92878935e-19],
[-1.50493730e-19, -4.46122155e-37, -1.55150365e-18,
9.98736570e-02, 0.00000000e+00, -8.22950122e-18,
-1.50493730e-19, 1.32252913e-33, 8.87128656e-18,
-6.42929783e-20, 1.92878935e-19, 1.72794621e-01]])
Or alternatively the hesse matrix:
mat = get_hesse_matrix(phonopy=phonopy_obj)
mat
array([[ 4.50127147e-02, -4.71227770e-32, 3.81802312e-32, ...,
-6.63514215e-05, 8.82979632e-06, 5.93920136e-05],
[-0.00000000e+00, 4.50127147e-02, -0.00000000e+00, ...,
8.82979632e-06, -6.63514215e-05, 5.93920136e-05],
[-0.00000000e+00, -0.00000000e+00, 4.50127147e-02, ...,
5.93659141e-05, 5.93659141e-05, 1.73512126e-05],
...,
[-6.63514215e-05, 8.82979632e-06, 5.93920136e-05, ...,
4.50127147e-02, 0.00000000e+00, 0.00000000e+00],
[ 8.82979632e-06, -6.63514215e-05, 5.93920136e-05, ...,
3.81802312e-32, 4.50127147e-02, -4.71227770e-32],
[ 5.93659141e-05, 5.93659141e-05, 1.73512126e-05, ...,
0.00000000e+00, 0.00000000e+00, 4.50127147e-02]])
Finally, also the function to calculate the band structure is directly available on the PhonopyWorkflow:
band_structure = get_band_structure(
phonopy=phonopy_obj,
npoints=101,
with_eigenvectors=False,
with_group_velocities=False
)
This band structure can also be visualised using the built-in plotting function:
plot_band_structure(
phonopy=phonopy_obj,
)
<Axes: title={'center': 'Bandstructure'}, xlabel='Bandpath', ylabel='Frequency [THz]'>
Just like the desnsity of states which can be plotted using:
plot_dos(phonopy_dict=phonopy_dict)
<Axes: title={'center': 'Phonon DOS vs Energy'}, xlabel='Frequency [THz]', ylabel='DOS'>
Quasi-harmonic Approximation#
To include the volume expansion with finite temperature the atomistics package implements the QuasiHarmonicWorkflow:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import (
get_tasks_for_quasi_harmonic_approximation,
analyse_results_for_quasi_harmonic_approximation,
get_thermal_properties_for_quasi_harmonic_approximation,
)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
task_dict, qh_dict = get_tasks_for_quasi_harmonic_approximation(
structure=bulk("Al", cubic=True),
num_points=11,
vol_range=0.05,
interaction_range=10,
displacement=0.01,
)
result_dict = evaluate_with_lammpslib(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
eng_internal_dict, phonopy_collect_dict = analyse_results_for_quasi_harmonic_approximation(
qh_dict=qh_dict,
output_dict=result_dict,
dos_mesh=20,
number_of_snapshots=None,
)
The QuasiHarmonicWorkflow is a combination of the EnergyVolumeCurveWorkflow and the PhonopyWorkflow. Consequently,
the inputs are a superset of the inputs of these two workflows.
Based on the QuasiHarmonicWorkflow the thermal properties can be calculated:
tp_dict = get_thermal_properties_for_quasi_harmonic_approximation(
eng_internal_dict=eng_internal_dict,
qh_dict=qh_dict,
task_dict=task_dict,
fit_type="polynomial",
fit_order=3,
t_min=1,
t_max=1500,
t_step=50,
temperatures=None,
cutoff_frequency=None,
pretend_real=False,
band_indices=None,
is_projection=False,
quantum_mechanical=True,
)
print(tp_dict)
{'temperatures': array([1.000e+00, 5.100e+01, 1.010e+02, 1.510e+02, 2.010e+02, 2.510e+02,
3.010e+02, 3.510e+02, 4.010e+02, 4.510e+02, 5.010e+02, 5.510e+02,
6.010e+02, 6.510e+02, 7.010e+02, 7.510e+02, 8.010e+02, 8.510e+02,
9.010e+02, 9.510e+02, 1.001e+03, 1.051e+03, 1.101e+03, 1.151e+03,
1.201e+03, 1.251e+03, 1.301e+03, 1.351e+03, 1.401e+03, 1.451e+03,
1.501e+03]), 'volumes': [66.71710763927877, 66.72177216699532, 66.75880304566003, 66.82252047263998, 66.89849494959407, 66.9796340113938, 67.06260790999805, 67.14571406293553, 67.2280041257585, 67.30891433770049, 67.38809533263728, 67.46532691224263, 67.54047194450703, 67.61344961443127, 67.68421888050432, 67.75276763760674, 67.8191052501301, 67.88325718225146, 67.94526099997334, 68.00516331350357, 68.06301739228013, 68.11888127927476, 68.17281628743957, 68.22488579579766, 68.27515428480672, 68.32368656530373, 68.3705471654409, 68.41579984728268, 68.45950723010068, 68.50173050155405, 68.54252920118518], 'free_energy': array([ 0.14903662, 0.14826796, 0.13934608, 0.1169922 , 0.08193524,
0.03597463, -0.01929655, -0.08261538, -0.15299036, -0.22963345,
-0.31190757, -0.39928889, -0.49134004, -0.5876908 , -0.68802399,
-0.79206502, -0.89957393, -1.01033927, -1.12417341, -1.24090869,
-1.36039449, -1.48249475, -1.60708597, -1.73405557, -1.86330055,
-1.99472628, -2.12824554, -2.26377775, -2.40124816, -2.54058732,
-2.6817305 ]), 'entropy': array([1.02970669e-08, 5.98072651e+00, 2.96865053e+01, 5.55852668e+01,
7.82409727e+01, 9.75218995e+01, 1.14065625e+02, 1.28465273e+02,
1.41174728e+02, 1.52530436e+02, 1.62783056e+02, 1.72122199e+02,
1.80693841e+02, 1.88612312e+02, 1.95968600e+02, 2.02836175e+02,
2.09275150e+02, 2.15335286e+02, 2.21058222e+02, 2.26479132e+02,
2.31627991e+02, 2.36530543e+02, 2.41209060e+02, 2.45682935e+02,
2.49969156e+02, 2.54082691e+02, 2.58036788e+02, 2.61843234e+02,
2.65512561e+02, 2.69054215e+02, 2.72476702e+02]), 'heat_capacity': array([1.67065980e-07, 1.73540235e+01, 5.38037700e+01, 7.36871465e+01,
8.35644372e+01, 8.88841670e+01, 9.20085315e+01, 9.39792227e+01,
9.52946945e+01, 9.62133951e+01, 9.68788951e+01, 9.73756862e+01,
9.77559504e+01, 9.80532534e+01, 9.82899463e+01, 9.84813619e+01,
9.86382931e+01, 9.87685101e+01, 9.88777197e+01, 9.89701872e+01,
9.90491516e+01, 9.91171073e+01, 9.91760000e+01, 9.92273651e+01,
9.92724273e+01, 9.93121724e+01, 9.93474017e+01, 9.93787712e+01,
9.94068224e+01, 9.94320054e+01, 9.94546967e+01])}
This requires the same inputs as the calculation of the thermal properties get_thermal_properties() with the
PhonopyWorkflow. The additional parameter quantum_mechanical specifies whether the classical harmonic oscillator or
the quantum mechanical harmonic oscillator is used to calculate the free energy.
And finally also the thermal expansion can be calculated:
tp_dict = get_thermal_properties_for_quasi_harmonic_approximation(
eng_internal_dict=eng_internal_dict,
qh_dict=qh_dict,
task_dict=task_dict,
fit_type="polynomial",
fit_order=3,
t_min=1,
t_max=1500,
t_step=50,
temperatures=None,
cutoff_frequency=None,
pretend_real=False,
band_indices=None,
is_projection=False,
quantum_mechanical=True,
output_keys=["free_energy", "temperatures", "volumes"],
)
temperatures, volumes = tp_dict["temperatures"], tp_dict["volumes"]
Structure Optimization#
In analogy to the molecular dynamics calculation also the structure optimization could in principle be defined inside
the simulation code or on the python level. Still currently the atomistics package only supports the structure
optimization defined inside the simulation codes.
Volume and Positions#
To optimize both the volume of the supercell as well as the positions inside the supercell the atomistics package
implements the optimize_positions_and_volume() workflow:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import optimize_positions_and_volume
structure = bulk("Al", a=4.0, cubic=True)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = evaluate_with_lammpslib(
task_dict=optimize_positions_and_volume(structure=structure),
potential_dataframe=potential_dataframe,
)
structure_opt = result_dict["structure_with_optimized_positions_and_volume"]
structure_opt
Atoms(symbols='Al4', pbc=True, cell=[4.05000466219724, 4.05000466219724, 4.05000466219724])
The result is the optimized atomistic structure as part of the result dictionary.
Positions#
The optimization of the positions inside the supercell without the optimization of the supercell volume is possible with
the optimize_positions() workflow:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammpslib, get_potential_by_name
from atomistics.workflows import optimize_positions
structure = bulk("Al", a=4.0, cubic=True)
potential_dataframe = get_potential_by_name(
potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", resource_path="static/lammps"
)
result_dict = evaluate_with_lammpslib(
task_dict=optimize_positions(structure=structure),
potential_dataframe=potential_dataframe,
)
structure_opt = result_dict["structure_with_optimized_positions"]
structure_opt
Atoms(symbols='Al4', pbc=True, cell=[4.0, 4.0, 4.0])
The result is the optimized atomistic structure as part of the result dictionary.