In [1]:
from atomistics.calculators import evaluate_with_lammps, get_potential_by_name
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
result_dict = evaluate_with_lammps(
task_dict={},
potential_dataframe=potential_dataframe,
)
The interatomic potential for Aluminium from Mishin named 1999--Mishin-Y--Al--LAMMPS--ipr1 is used in the evaluation
with LAMMPS evaluate_with_lammps().
Elastic Matrix¶
The elastic constants and elastic moduli can be calculated using the ElasticMatrixWorkflow:
In [2]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, get_potential_by_name
from atomistics.workflows import ElasticMatrixWorkflow
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
workflow = ElasticMatrixWorkflow(
structure=bulk("Al", cubic=True),
num_of_point=5,
eps_range=0.005,
sqrt_eta=True,
fit_order=2,
)
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
fit_dict = workflow.analyse_structures(output_dict=result_dict)
print(fit_dict)
OrderedDict([('SGN', 225), ('v0', 66.43012500000002), ('LC', 'CI'), ('Lag_strain_list', ['01', '08', '23']), ('epss', array([-0.005 , -0.0025, 0. , 0.0025, 0.005 ])), ('e0', -13.439999952539933), ('strain_energy', [[(-0.005, -13.436318571477456), (-0.0025, -13.439078837682322), (0.0, -13.439999952539933), (0.0024999999999999996, -13.439085818358748), (0.005, -13.436366000141689)], [(-0.005, -13.438173590018806), (-0.0025, -13.43954407461901), (0.0, -13.439999952539933), (0.0024999999999999996, -13.439548790708864), (0.005, -13.438205313732144)], [(-0.005, -13.436741916886653), (-0.0025, -13.439195456175232), (0.0, -13.439999952539933), (0.0024999999999999996, -13.439213481749942), (0.005, -13.436885676373926)]]), ('C', array([[114.10393023, 60.51098897, 60.51098897, 0. ,
0. , 0. ],
[ 60.51098897, 114.10393023, 60.51098897, 0. ,
0. , 0. ],
[ 60.51098897, 60.51098897, 114.10393023, 0. ,
0. , 0. ],
[ 0. , 0. , 0. , 51.23931149,
0. , 0. ],
[ 0. , 0. , 0. , 0. ,
51.23931149, 0. ],
[ 0. , 0. , 0. , 0. ,
0. , 51.23931149]])), ('A2', array([2.20131073, 1.0898606 , 1.91886377])), ('BV', 78.37530272473929), ('GV', 41.462175146677424), ('EV', 105.74025607889799), ('nuV', 0.2751412064710683), ('S', 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.01951627, 0. ,
0. ],
[ 0. , 0. , 0. , 0. , 0.01951627,
0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0.01951627]])), ('BR', 78.37530272473931), ('GR', 37.54162684596518), ('ER', 97.1183728107761), ('nuR', 0.29347581564934205), ('BH', 78.3753027247393), ('GH', 39.501900996321304), ('EH', 101.46008564559224), ('nuH', 0.2842430754793411), ('AVR', 4.962480541224269), ('C_eigval', EigResult(eigenvalues=array([ 53.59294126, 235.12590817, 53.59294126, 51.23931149,
51.23931149, 51.23931149]), eigenvectors=array([[-0.81649658, 0.57735027, 0.11541902, 0. , 0. ,
0. ],
[ 0.40824829, 0.57735027, -0.75771582, 0. , 0. ,
0. ],
[ 0.40824829, 0.57735027, 0.6422968 , 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.
In [3]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, get_potential_by_name
from atomistics.workflows import EnergyVolumeCurveWorkflow
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
workflow = EnergyVolumeCurveWorkflow(
structure=bulk("Al", cubic=True),
num_points=11,
fit_type="polynomial",
fit_order=3,
vol_range=0.05,
axes=("x", "y", "z"),
strains=None,
)
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
fit_dict = workflow.analyse_structures(output_dict=result_dict)
print(fit_dict)
{'poly_fit': array([-4.17645808e-05, 1.19746500e-02, -1.03803906e+00, 1.49168639e+01]), 'fit_type': 'polynomial', 'fit_order': 3, 'volume_eq': 66.43019853103964, 'energy_eq': -13.43996804374383, 'bulkmodul_eq': 77.7250135953191, 'b_prime_eq': 1.279502459079921, 'least_square_error': 3.225313797039607e-10, 'volume': [63.10861874999998, 63.77291999999998, 64.43722124999998, 65.1015225, 65.76582375000004, 66.43012500000002, 67.09442624999994, 67.75872750000002, 68.42302874999999, 69.08732999999997, 69.75163125000002], 'energy': [-13.398169481534445, -13.413389552957456, -13.425112589013958, -13.433411420804067, -13.438357630783006, -13.439999952539933, -13.438383476946305, -13.433607982916406, -13.425774537190858, -13.414961805921427, -13.401233093668836]}
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:
In [4]:
tp_dict = workflow.get_thermal_properties(
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]), '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.08219062, 18.62461552, 38.05446426,
57.6693229 , 75.37710506, 90.99476554, 104.78762778,
117.06473011, 128.09164494, 138.08127289, 147.20167195,
155.58579193, 163.33970927, 170.54896552, 177.28330938,
183.60022562, 189.54757244, 195.16556897, 200.48830826,
205.54492122, 210.36048158, 214.95671661, 219.35257076,
223.56465688, 227.60762034, 231.49443548, 235.23664867,
238.84457908, 242.32748555, 244.0403182 ]), 'heat_capacity': array([8.65067172e-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]), '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])}
Or alternatively directly calculate the thermal expansion:
In [5]:
temperatures, volumes = workflow.get_thermal_expansion(
output_dict=result_dict,
t_min=1,
t_max=1500,
t_step=50,
temperatures=None,
)
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.
In [6]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_langevin_with_lammps,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
result_dict = calc_molecular_dynamics_langevin_with_lammps(
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=("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:
In [7]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_nvt_with_lammps,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
result_dict = calc_molecular_dynamics_nvt_with_lammps(
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=("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:
In [8]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_npt_with_lammps,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
result_dict = calc_molecular_dynamics_npt_with_lammps(
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=("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:
In [9]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_nph_with_lammps,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
result_dict = calc_molecular_dynamics_nph_with_lammps(
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=("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:
In [10]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_thermal_expansion_with_lammps,
evaluate_with_lammps,
get_potential_by_name,
)
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
temperatures_md, volumes_md = calc_molecular_dynamics_thermal_expansion_with_lammps(
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,
)
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:
In [11]:
from ase.build import bulk
from atomistics.calculators import (
calc_molecular_dynamics_phonons_with_lammps,
evaluate_with_lammps,
)
from atomistics.workflows import optimize_positions_and_volume, PhonopyWorkflow
from dynaphopy import Quasiparticle
import pandas
from phonopy.units import VaspToTHz
import spglib
structure_bulk = bulk("Si", cubic=True)
potential_dataframe = get_potential_by_name(
potential_name='1988--Tersoff-J--Si-c--LAMMPS--ipr1'
)
The first step is optimizing the Silicon diamond structure to match the lattice specifications implemented in the Tersoff interatomic potential:
In [12]:
task_dict = optimize_positions_and_volume(structure=structure_bulk)
result_dict = evaluate_with_lammps(
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.
In [13]:
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)
workflow = PhonopyWorkflow(
structure=structure_ase,
interaction_range=10,
factor=VaspToTHz,
displacement=0.01,
dos_mesh=20,
primitive_matrix=primitive_matrix,
number_of_snapshots=None,
)
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
workflow.analyse_structures(output_dict=result_dict)
Out[13]:
({'qpoints': array([[0.025, 0.025, 0.025],
[0.075, 0.025, 0.025],
[0.125, 0.025, 0.025],
...,
[0.525, 0.525, 0.425],
[0.475, 0.475, 0.475],
[0.525, 0.475, 0.475]]),
'weights': array([ 2, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12,
12, 12, 6, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 6, 12,
12, 12, 12, 6, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12,
12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12,
12, 6, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 6, 12, 12, 12, 6,
12, 12, 6, 12, 6, 2, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12,
12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 6, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12,
12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12,
12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12,
6, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 12,
12, 6, 12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 6, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 6,
12, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12, 12, 6, 12, 12,
12, 12, 12, 6, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6, 6,
6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 12, 12,
6, 12, 12, 12, 12, 6, 12, 12, 6, 12, 6, 2, 6, 6, 6, 6, 6,
6, 12, 12, 12, 12, 6, 12, 12, 12, 6, 12, 6, 2, 6, 6, 6, 6,
12, 12, 6, 2, 6]),
'frequencies': array([[ 0.35167322, 0.35167322, 0.71941397, 16.05999349, 16.06435417,
16.06435417],
[ 0.7810313 , 1.04550359, 1.76442286, 16.01908158, 16.04238401,
16.0474694 ],
[ 1.43010566, 1.69629248, 3.13109956, 15.91248433, 15.99513531,
16.00361669],
...,
[ 4.89126604, 5.30105813, 11.18409265, 13.02666668, 15.38025565,
15.46187631],
[ 4.65215064, 4.65215064, 11.21184039, 13.22967807, 15.43087981,
15.43087981],
[ 4.71432047, 4.84620075, 11.26998698, 13.12729988, 15.41640899,
15.43613968]]),
'eigenvectors': None,
'group_velocities': None},
{'frequency_points': array([-1.21959488e+00, -1.12531879e+00, -1.03104271e+00, -9.36766620e-01,
-8.42490534e-01, -7.48214449e-01, -6.53938363e-01, -5.59662277e-01,
-4.65386192e-01, -3.71110106e-01, -2.76834020e-01, -1.82557935e-01,
-8.82818488e-02, 5.99423689e-03, 1.00270323e-01, 1.94546408e-01,
2.88822494e-01, 3.83098580e-01, 4.77374665e-01, 5.71650751e-01,
6.65926837e-01, 7.60202922e-01, 8.54479008e-01, 9.48755094e-01,
1.04303118e+00, 1.13730727e+00, 1.23158335e+00, 1.32585944e+00,
1.42013552e+00, 1.51441161e+00, 1.60868769e+00, 1.70296378e+00,
1.79723987e+00, 1.89151595e+00, 1.98579204e+00, 2.08006812e+00,
2.17434421e+00, 2.26862029e+00, 2.36289638e+00, 2.45717246e+00,
2.55144855e+00, 2.64572464e+00, 2.74000072e+00, 2.83427681e+00,
2.92855289e+00, 3.02282898e+00, 3.11710506e+00, 3.21138115e+00,
3.30565724e+00, 3.39993332e+00, 3.49420941e+00, 3.58848549e+00,
3.68276158e+00, 3.77703766e+00, 3.87131375e+00, 3.96558984e+00,
4.05986592e+00, 4.15414201e+00, 4.24841809e+00, 4.34269418e+00,
4.43697026e+00, 4.53124635e+00, 4.62552244e+00, 4.71979852e+00,
4.81407461e+00, 4.90835069e+00, 5.00262678e+00, 5.09690286e+00,
5.19117895e+00, 5.28545504e+00, 5.37973112e+00, 5.47400721e+00,
5.56828329e+00, 5.66255938e+00, 5.75683546e+00, 5.85111155e+00,
5.94538764e+00, 6.03966372e+00, 6.13393981e+00, 6.22821589e+00,
6.32249198e+00, 6.41676806e+00, 6.51104415e+00, 6.60532024e+00,
6.69959632e+00, 6.79387241e+00, 6.88814849e+00, 6.98242458e+00,
7.07670066e+00, 7.17097675e+00, 7.26525284e+00, 7.35952892e+00,
7.45380501e+00, 7.54808109e+00, 7.64235718e+00, 7.73663326e+00,
7.83090935e+00, 7.92518544e+00, 8.01946152e+00, 8.11373761e+00,
8.20801369e+00, 8.30228978e+00, 8.39656586e+00, 8.49084195e+00,
8.58511804e+00, 8.67939412e+00, 8.77367021e+00, 8.86794629e+00,
8.96222238e+00, 9.05649846e+00, 9.15077455e+00, 9.24505064e+00,
9.33932672e+00, 9.43360281e+00, 9.52787889e+00, 9.62215498e+00,
9.71643106e+00, 9.81070715e+00, 9.90498323e+00, 9.99925932e+00,
1.00935354e+01, 1.01878115e+01, 1.02820876e+01, 1.03763637e+01,
1.04706397e+01, 1.05649158e+01, 1.06591919e+01, 1.07534680e+01,
1.08477441e+01, 1.09420202e+01, 1.10362963e+01, 1.11305723e+01,
1.12248484e+01, 1.13191245e+01, 1.14134006e+01, 1.15076767e+01,
1.16019528e+01, 1.16962289e+01, 1.17905049e+01, 1.18847810e+01,
1.19790571e+01, 1.20733332e+01, 1.21676093e+01, 1.22618854e+01,
1.23561615e+01, 1.24504375e+01, 1.25447136e+01, 1.26389897e+01,
1.27332658e+01, 1.28275419e+01, 1.29218180e+01, 1.30160941e+01,
1.31103701e+01, 1.32046462e+01, 1.32989223e+01, 1.33931984e+01,
1.34874745e+01, 1.35817506e+01, 1.36760267e+01, 1.37703027e+01,
1.38645788e+01, 1.39588549e+01, 1.40531310e+01, 1.41474071e+01,
1.42416832e+01, 1.43359593e+01, 1.44302353e+01, 1.45245114e+01,
1.46187875e+01, 1.47130636e+01, 1.48073397e+01, 1.49016158e+01,
1.49958919e+01, 1.50901679e+01, 1.51844440e+01, 1.52787201e+01,
1.53729962e+01, 1.54672723e+01, 1.55615484e+01, 1.56558245e+01,
1.57501005e+01, 1.58443766e+01, 1.59386527e+01, 1.60329288e+01,
1.61272049e+01, 1.62214810e+01, 1.63157571e+01, 1.64100331e+01,
1.65043092e+01, 1.65985853e+01, 1.66928614e+01, 1.67871375e+01,
1.68814136e+01, 1.69756897e+01, 1.70699657e+01, 1.71642418e+01,
1.72585179e+01, 1.73527940e+01, 1.74470701e+01, 1.75413462e+01,
1.76356223e+01]),
'total_dos': array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 3.17526949e-04, 1.35171618e-03, 2.50831797e-03,
3.78733233e-03, 5.26373813e-03, 6.90328284e-03, 8.78451412e-03,
1.08845716e-02, 1.31202482e-02, 1.58239881e-02, 1.84817583e-02,
2.14099537e-02, 2.44295300e-02, 2.78000795e-02, 3.14939706e-02,
3.52413320e-02, 3.92909859e-02, 4.36423577e-02, 4.83495504e-02,
5.31903107e-02, 5.81462801e-02, 6.35669356e-02, 6.94844360e-02,
7.55171874e-02, 8.19720137e-02, 8.86760193e-02, 9.58023138e-02,
1.03099606e-01, 1.11366958e-01, 1.19969573e-01, 1.28656383e-01,
1.38210214e-01, 1.48087177e-01, 1.58993039e-01, 1.70577194e-01,
1.82773665e-01, 1.96102840e-01, 2.10318186e-01, 2.25609695e-01,
2.42775583e-01, 2.61718244e-01, 2.82966510e-01, 3.07071183e-01,
3.35572804e-01, 3.71503640e-01, 4.23538191e-01, 5.01245367e-01,
5.06757340e-01, 5.10673514e-01, 5.13884731e-01, 5.18635687e-01,
5.22215630e-01, 5.26790868e-01, 5.30317993e-01, 5.33703746e-01,
5.38149099e-01, 5.41609848e-01, 5.45017573e-01, 5.48177732e-01,
5.52308667e-01, 5.54935999e-01, 5.58609041e-01, 5.61624143e-01,
5.64469403e-01, 5.68022742e-01, 5.70890330e-01, 5.72980739e-01,
5.77439994e-01, 5.83268678e-01, 5.49870464e-01, 4.81961887e-01,
4.49243839e-01, 4.26254220e-01, 4.10187378e-01, 3.95199652e-01,
3.98624284e-01, 4.05995249e-01, 4.25103229e-01, 4.54252120e-01,
4.65563502e-01, 3.65773926e-01, 2.85541037e-01, 2.23526600e-01,
1.44473263e-01, 1.04439534e-01, 1.08917911e-01, 1.13417776e-01,
1.18953226e-01, 1.24441198e-01, 1.29621905e-01, 1.36164980e-01,
1.42605861e-01, 1.49045426e-01, 1.56851391e-01, 1.64582559e-01,
1.72682584e-01, 1.82141611e-01, 1.91321788e-01, 2.02261449e-01,
2.13762731e-01, 2.26135950e-01, 2.40393404e-01, 2.55442497e-01,
2.74046636e-01, 2.93902607e-01, 3.19085171e-01, 3.48960866e-01,
3.91473391e-01, 4.54354038e-01, 5.99226889e-01, 6.56603266e-01,
4.43486451e-01, 3.58148235e-01, 2.92612110e-01, 2.34458733e-01,
1.57664249e-01, 9.44996839e-02, 7.64989495e-02, 6.21782198e-02,
7.46273535e-02, 8.86461614e-02, 1.06333999e-01, 1.36147948e-01,
1.87831258e-01, 2.67221370e-01, 2.89615750e-01, 3.09688343e-01,
3.30425521e-01, 3.30668402e-01, 3.40062237e-01, 3.45284499e-01,
3.55917538e-01, 3.60077500e-01, 3.69186627e-01, 3.71502854e-01,
3.72120651e-01, 3.48181936e-01, 3.16317708e-01, 2.98435672e-01,
2.84862196e-01, 2.74096441e-01, 2.64721157e-01, 2.56360104e-01,
2.48761334e-01, 2.41852996e-01, 2.35018844e-01, 2.28859604e-01,
2.22571780e-01, 2.16655364e-01, 2.10807826e-01, 2.04860365e-01,
1.99221497e-01, 1.93100666e-01, 1.87355107e-01, 2.53244401e-01,
8.43177422e-01, 1.43196664e+00, 3.18170027e+00, 2.79552678e+00,
3.13543115e+00, 4.55809708e+00, 2.72396129e+00, 1.52338820e+00,
1.06733102e+00, 7.67955185e-01, 5.18138351e-01, 2.28096624e-01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00])})
The calcualtion of the finite temperature phonons starts by computing the molecular dynamics trajectory using the
calc_molecular_dynamics_phonons_with_lammps() function. This function is internally linked to DynaPhoPy
to return an dynaphopy.dynamics.Dynamics object:
In [14]:
trajectory = calc_molecular_dynamics_phonons_with_lammps(
structure_ase=structure_ase,
potential_dataframe=potential_dataframe,
force_constants=workflow.phonopy.get_force_constants(),
phonopy_unitcell=workflow.phonopy.get_unitcell(),
phonopy_primitive_matrix=workflow.phonopy.get_primitive_matrix(),
phonopy_supercell_matrix=workflow.phonopy.get_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:
In [15]:
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): [2.18910938e-06 2.19854601e-06 2.20383682e-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.503412 THz Position 0.030228 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 779538375233.310425 Frequency shift 0.030226 THz Peak # 2 ---------------------------------------------- Width 0.503412 THz Position 0.030228 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 779538375233.310425 Frequency shift 0.030226 THz Peak # 3 ---------------------------------------------- Width 0.503412 THz Position 0.030228 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 779538375233.310425 Frequency shift 0.030226 THz Peak # 4 ---------------------------------------------- Width 0.786715 THz Position 15.561772 THz Area (<K>) (Lorentzian) 0.014497 eV Area (<K>) (Total) 0.013722 eV <|dQ/dt|^2> 0.028993 eV Occupation number 2.330539 Fit temperature 332.921389 K Base line -0.000016 eV * ps Maximum height 0.011731 eV * ps Fitting global error 0.033291 Frequency shift -0.506127 THz Peak # 5 ---------------------------------------------- Width 0.786715 THz Position 15.561772 THz Area (<K>) (Lorentzian) 0.014497 eV Area (<K>) (Total) 0.013722 eV <|dQ/dt|^2> 0.028993 eV Occupation number 2.330539 Fit temperature 332.921389 K Base line -0.000016 eV * ps Maximum height 0.011731 eV * ps Fitting global error 0.033291 Frequency shift -0.506127 THz Peak # 6 ---------------------------------------------- Width 0.786715 THz Position 15.561772 THz Area (<K>) (Lorentzian) 0.014497 eV Area (<K>) (Total) 0.013722 eV <|dQ/dt|^2> 0.028993 eV Occupation number 2.330539 Fit temperature 332.921389 K Base line -0.000016 eV * ps Maximum height 0.011731 eV * ps Fitting global error 0.033291 Frequency shift -0.506127 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.520811 THz Position 4.512520 THz Area (<K>) (Lorentzian) 0.018115 eV Area (<K>) (Total) 0.016752 eV <|dQ/dt|^2> 0.036230 eV Occupation number 11.697748 Fit temperature 420.192477 K Base line -0.000044 eV * ps Maximum height 0.022143 eV * ps Fitting global error 0.016927 Frequency shift -0.151454 THz Peak # 2 ---------------------------------------------- Width 0.520811 THz Position 4.512520 THz Area (<K>) (Lorentzian) 0.018115 eV Area (<K>) (Total) 0.016752 eV <|dQ/dt|^2> 0.036230 eV Occupation number 11.697748 Fit temperature 420.192477 K Base line -0.000044 eV * ps Maximum height 0.022143 eV * ps Fitting global error 0.016927 Frequency shift -0.151454 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.575794 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.942499 THz Position 15.067970 THz Area (<K>) (Lorentzian) 0.022568 eV Area (<K>) (Total) 0.024331 eV <|dQ/dt|^2> 0.045136 eV Occupation number 4.050999 Fit temperature 521.672324 K Base line 0.000120 eV * ps Maximum height 0.015244 eV * ps Fitting global error 0.020341 Frequency shift -0.487709 THz Peak # 6 ---------------------------------------------- Width 0.942499 THz Position 15.067970 THz Area (<K>) (Lorentzian) 0.022568 eV Area (<K>) (Total) 0.024331 eV <|dQ/dt|^2> 0.045136 eV Occupation number 4.050999 Fit temperature 521.672324 K Base line 0.000120 eV * ps Maximum height 0.015244 eV * ps Fitting global error 0.020341 Frequency shift -0.487709 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.579259 THz Position 6.561490 THz Area (<K>) (Lorentzian) 0.025327 eV Area (<K>) (Total) 0.027369 eV <|dQ/dt|^2> 0.050654 eV Occupation number 11.228659 Fit temperature 587.462943 K Base line 0.000121 eV * ps Maximum height 0.027835 eV * ps Fitting global error 0.025639 Frequency shift -0.333845 THz Peak # 2 ---------------------------------------------- Width 0.579259 THz Position 6.561490 THz Area (<K>) (Lorentzian) 0.025327 eV Area (<K>) (Total) 0.027369 eV <|dQ/dt|^2> 0.050654 eV Occupation number 11.228659 Fit temperature 587.462943 K Base line 0.000121 eV * ps Maximum height 0.027835 eV * ps Fitting global error 0.025639 Frequency shift -0.333845 THz Peak # 3 ---------------------------------------------- Width 0.605260 THz Position 11.933579 THz Area (<K>) (Lorentzian) 0.030516 eV Area (<K>) (Total) 0.034730 eV <|dQ/dt|^2> 0.061032 eV Occupation number 7.270035 Fit temperature 707.271379 K Base line 0.000225 eV * ps Maximum height 0.032097 eV * ps Fitting global error 0.023688 Frequency shift -0.258211 THz Peak # 4 ---------------------------------------------- Width 0.605260 THz Position 11.933579 THz Area (<K>) (Lorentzian) 0.030516 eV Area (<K>) (Total) 0.034730 eV <|dQ/dt|^2> 0.061032 eV Occupation number 7.270035 Fit temperature 707.271379 K Base line 0.000225 eV * ps Maximum height 0.032097 eV * ps Fitting global error 0.023688 Frequency shift -0.258211 THz Peak # 5 ---------------------------------------------- Width 0.634083 THz Position 14.446261 THz Area (<K>) (Lorentzian) 0.042334 eV Area (<K>) (Total) 0.048339 eV <|dQ/dt|^2> 0.084669 eV Occupation number 8.404323 Fit temperature 981.504235 K Base line 0.000327 eV * ps Maximum height 0.042504 eV * ps Fitting global error 0.015367 Frequency shift -0.444694 THz Peak # 6 ---------------------------------------------- Width 0.634083 THz Position 14.446261 THz Area (<K>) (Lorentzian) 0.042334 eV Area (<K>) (Total) 0.048339 eV <|dQ/dt|^2> 0.084669 eV Occupation number 8.404323 Fit temperature 981.504235 K Base line 0.000327 eV * ps Maximum height 0.042504 eV * ps Fitting global error 0.015367 Frequency shift -0.444694 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.513828 THz Position 4.501741 THz Area (<K>) (Lorentzian) 0.024073 eV Area (<K>) (Total) 0.021184 eV <|dQ/dt|^2> 0.048147 eV Occupation number 15.748758 Fit temperature 558.542799 K Base line -0.000110 eV * ps Maximum height 0.029826 eV * ps Fitting global error 0.014477 Frequency shift -0.162232 THz Peak # 2 ---------------------------------------------- Width 0.513828 THz Position 4.501741 THz Area (<K>) (Lorentzian) 0.024073 eV Area (<K>) (Total) 0.021184 eV <|dQ/dt|^2> 0.048147 eV Occupation number 15.748758 Fit temperature 558.542799 K Base line -0.000110 eV * ps Maximum height 0.029826 eV * ps Fitting global error 0.014477 Frequency shift -0.162232 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.018837 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.439740 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.693495 THz Position 15.046973 THz Area (<K>) (Lorentzian) 0.030097 eV Area (<K>) (Total) 0.029819 eV <|dQ/dt|^2> 0.060195 eV Occupation number 5.577759 Fit temperature 696.952048 K Base line 0.000023 eV * ps Maximum height 0.027629 eV * ps Fitting global error 0.016708 Frequency shift -0.508706 THz Peak # 6 ---------------------------------------------- Width 0.693495 THz Position 15.046973 THz Area (<K>) (Lorentzian) 0.030097 eV Area (<K>) (Total) 0.029819 eV <|dQ/dt|^2> 0.060195 eV Occupation number 5.577759 Fit temperature 696.952048 K Base line 0.000023 eV * ps Maximum height 0.027629 eV * ps Fitting global error 0.016708 Frequency shift -0.508706 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... 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.528069 THz Position 4.477183 THz Area (<K>) (Lorentzian) 0.042227 eV Area (<K>) (Total) 0.040708 eV <|dQ/dt|^2> 0.084453 eV Occupation number 28.158071 Fit temperature 979.942726 K Base line -0.000024 eV * ps Maximum height 0.050907 eV * ps Fitting global error 0.012211 Frequency shift -0.190696 THz Peak # 2 ---------------------------------------------- Width 0.890764 THz Position 6.695151 THz Area (<K>) (Lorentzian) 0.080890 eV Area (<K>) (Total) 0.093342 eV <|dQ/dt|^2> 0.161780 eV Occupation number 36.211198 Fit temperature 1877.262476 K Base line 0.000706 eV * ps Maximum height 0.057811 eV * ps Fitting global error 0.020823 Frequency shift -0.265939 THz Peak # 3 ---------------------------------------------- Width 0.879866 THz Position 8.803246 THz Area (<K>) (Lorentzian) 0.043886 eV Area (<K>) (Total) 0.052217 eV <|dQ/dt|^2> 0.087772 eV Occupation number 14.647678 Fit temperature 1018.178750 K Base line 0.000449 eV * ps Maximum height 0.031753 eV * ps Fitting global error 0.028094 Frequency shift -0.202601 THz Peak # 4 ---------------------------------------------- Width 0.757650 THz Position 13.371395 THz Area (<K>) (Lorentzian) 0.021906 eV Area (<K>) (Total) 0.026391 eV <|dQ/dt|^2> 0.043812 eV Occupation number 4.477924 Fit temperature 506.700042 K Base line 0.000237 eV * ps Maximum height 0.018407 eV * ps Fitting global error 0.037761 Frequency shift -0.353521 THz Peak # 5 ---------------------------------------------- Width 0.668122 THz Position 14.983805 THz Area (<K>) (Lorentzian) 0.023784 eV Area (<K>) (Total) 0.024310 eV <|dQ/dt|^2> 0.047568 eV Occupation number 4.323155 Fit temperature 550.026026 K Base line 0.000052 eV * ps Maximum height 0.022663 eV * ps Fitting global error 0.013903 Frequency shift -0.442641 THz Peak # 6 ---------------------------------------------- Width 0.793991 THz Position 15.147433 THz Area (<K>) (Lorentzian) 0.067893 eV Area (<K>) (Total) 0.078771 eV <|dQ/dt|^2> 0.135785 eV Occupation number 13.119087 Fit temperature 1575.015038 K Base line 0.000607 eV * ps Maximum height 0.054436 eV * ps Fitting global error 0.022076 Frequency shift -0.435323 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... 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.906814 THz Position 7.245967 THz Area (<K>) (Lorentzian) 0.086091 eV Area (<K>) (Total) 0.103019 eV <|dQ/dt|^2> 0.172181 eV Occupation number 35.601390 Fit temperature 1997.953650 K Base line 0.000922 eV * ps Maximum height 0.060439 eV * ps Fitting global error 0.020747 Frequency shift -0.296529 THz Peak # 2 ---------------------------------------------- Width 0.906814 THz Position 7.245967 THz Area (<K>) (Lorentzian) 0.086091 eV Area (<K>) (Total) 0.103019 eV <|dQ/dt|^2> 0.172181 eV Occupation number 35.601390 Fit temperature 1997.953650 K Base line 0.000922 eV * ps Maximum height 0.060439 eV * ps Fitting global error 0.020747 Frequency shift -0.296529 THz Peak # 3 ---------------------------------------------- Width 0.632257 THz Position 11.066196 THz Area (<K>) (Lorentzian) 0.021718 eV Area (<K>) (Total) 0.024990 eV <|dQ/dt|^2> 0.043436 eV Occupation number 5.463247 Fit temperature 502.867083 K Base line 0.000174 eV * ps Maximum height 0.021868 eV * ps Fitting global error 0.025864 Frequency shift -0.284125 THz Peak # 4 ---------------------------------------------- Width 0.632257 THz Position 11.066196 THz Area (<K>) (Lorentzian) 0.021718 eV Area (<K>) (Total) 0.024990 eV <|dQ/dt|^2> 0.043436 eV Occupation number 5.463247 Fit temperature 502.867083 K Base line 0.000174 eV * ps Maximum height 0.021868 eV * ps Fitting global error 0.025864 Frequency shift -0.284125 THz Peak # 5 ---------------------------------------------- Width 0.834964 THz Position 14.801288 THz Area (<K>) (Lorentzian) 0.039949 eV Area (<K>) (Total) 0.044302 eV <|dQ/dt|^2> 0.079898 eV Occupation number 7.701072 Fit temperature 926.027948 K Base line 0.000261 eV * ps Maximum height 0.030459 eV * ps Fitting global error 0.031087 Frequency shift -0.437050 THz Peak # 6 ---------------------------------------------- Width 0.834964 THz Position 14.801288 THz Area (<K>) (Lorentzian) 0.039949 eV Area (<K>) (Total) 0.044302 eV <|dQ/dt|^2> 0.079898 eV Occupation number 7.701072 Fit temperature 926.027948 K Base line 0.000261 eV * ps Maximum height 0.030459 eV * ps Fitting global error 0.031087 Frequency shift -0.437050 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... Harmonic frequencies (THz): [ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543] Peak # 1 ---------------------------------------------- Width 0.522556 THz Position 4.483775 THz Area (<K>) (Lorentzian) 0.037188 eV Area (<K>) (Total) 0.034635 eV <|dQ/dt|^2> 0.074376 eV Occupation number 24.701308 Fit temperature 862.984209 K Base line -0.000079 eV * ps Maximum height 0.045305 eV * ps Fitting global error 0.012393 Frequency shift -0.184104 THz Peak # 2 ---------------------------------------------- Width 0.898262 THz Position 6.707279 THz Area (<K>) (Lorentzian) 0.067887 eV Area (<K>) (Total) 0.079212 eV <|dQ/dt|^2> 0.135775 eV Occupation number 30.254428 Fit temperature 1575.464575 K Base line 0.000635 eV * ps Maximum height 0.048114 eV * ps Fitting global error 0.023773 Frequency shift -0.253812 THz Peak # 3 ---------------------------------------------- Width 0.835420 THz Position 8.803037 THz Area (<K>) (Lorentzian) 0.013150 eV Area (<K>) (Total) 0.016324 eV <|dQ/dt|^2> 0.026301 eV Occupation number 4.039095 Fit temperature 303.968668 K Base line 0.000166 eV * ps Maximum height 0.010021 eV * ps Fitting global error 0.059276 Frequency shift -0.202810 THz Peak # 4 ---------------------------------------------- Width 0.730369 THz Position 13.368389 THz Area (<K>) (Lorentzian) 0.016678 eV Area (<K>) (Total) 0.021238 eV <|dQ/dt|^2> 0.033357 eV Occupation number 3.290872 Fit temperature 384.834116 K Base line 0.000234 eV * ps Maximum height 0.014538 eV * ps Fitting global error 0.045590 Frequency shift -0.356527 THz Peak # 5 ---------------------------------------------- Width 0.674322 THz Position 14.951363 THz Area (<K>) (Lorentzian) 0.036155 eV Area (<K>) (Total) 0.038613 eV <|dQ/dt|^2> 0.072311 eV Occupation number 6.847775 Fit temperature 837.833759 K Base line 0.000157 eV * ps Maximum height 0.034134 eV * ps Fitting global error 0.014619 Frequency shift -0.475083 THz Peak # 6 ---------------------------------------------- Width 0.902049 THz Position 15.073779 THz Area (<K>) (Lorentzian) 0.021209 eV Area (<K>) (Total) 0.022490 eV <|dQ/dt|^2> 0.042418 eV Occupation number 3.775244 Fit temperature 489.986369 K Base line 0.000094 eV * ps Maximum height 0.014968 eV * ps Fitting global error 0.022499 Frequency shift -0.508977 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] 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... Peak # 1 ---------------------------------------------- Width 0.521991 THz Position 4.484158 THz Area (<K>) (Lorentzian) 0.036566 eV Area (<K>) (Total) 0.033926 eV <|dQ/dt|^2> 0.073131 eV Occupation number 24.277431 Fit temperature 848.537722 K Base line -0.000083 eV * ps Maximum height 0.044595 eV * ps Fitting global error 0.012530 Frequency shift -0.183721 THz Peak # 2 ---------------------------------------------- Width 0.899617 THz Position 6.704187 THz Area (<K>) (Lorentzian) 0.061915 eV Area (<K>) (Total) 0.071943 eV <|dQ/dt|^2> 0.123830 eV Occupation number 27.561609 Fit temperature 1436.830955 K Base line 0.000565 eV * ps Maximum height 0.043814 eV * ps Fitting global error 0.025259 Frequency shift -0.256903 THz Peak # 3 ---------------------------------------------- Width 0.840794 THz Position 8.777709 THz Area (<K>) (Lorentzian) 0.009473 eV Area (<K>) (Total) 0.011919 eV <|dQ/dt|^2> 0.018945 eV Occupation number 2.779076 Fit temperature 218.134958 K Base line 0.000127 eV * ps Maximum height 0.007172 eV * ps Fitting global error 0.078651 Frequency shift -0.228137 THz Peak # 4 ---------------------------------------------- Width 0.814357 THz Position 13.294667 THz Area (<K>) (Lorentzian) 0.012215 eV Area (<K>) (Total) 0.016240 eV <|dQ/dt|^2> 0.024429 eV Occupation number 2.291679 Fit temperature 280.431128 K Base line 0.000205 eV * ps Maximum height 0.009549 eV * ps Fitting global error 0.061800 Frequency shift -0.430249 THz Peak # 5 ---------------------------------------------- Width 0.706639 THz Position 14.920962 THz Area (<K>) (Lorentzian) 0.033379 eV Area (<K>) (Total) 0.034927 eV <|dQ/dt|^2> 0.066758 eV Occupation number 6.297384 Fit temperature 773.297147 K Base line 0.000113 eV * ps Maximum height 0.030072 eV * ps Fitting global error 0.021423 Frequency shift -0.505484 THz Peak # 6 ---------------------------------------------- Width 0.901748 THz Position 15.110122 THz Area (<K>) (Lorentzian) 0.021178 eV Area (<K>) (Total) 0.023049 eV <|dQ/dt|^2> 0.042356 eV Occupation number 3.758704 Fit temperature 489.249958 K Base line 0.000121 eV * ps Maximum height 0.014951 eV * ps Fitting global error 0.028888 Frequency shift -0.472634 THz 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] 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... Harmonic frequencies (THz): [ 4.66787904 6.96109048 9.00584683 13.72491589 15.42644585 15.58275543] Peak # 1 ---------------------------------------------- Width 0.521626 THz Position 4.484270 THz Area (<K>) (Lorentzian) 0.040369 eV Area (<K>) (Total) 0.037483 eV <|dQ/dt|^2> 0.080738 eV Occupation number 26.853838 Fit temperature 936.816652 K Base line -0.000091 eV * ps Maximum height 0.049268 eV * ps Fitting global error 0.011912 Frequency shift -0.183609 THz Peak # 2 ---------------------------------------------- Width 0.898998 THz Position 6.706088 THz Area (<K>) (Lorentzian) 0.060155 eV Area (<K>) (Total) 0.069989 eV <|dQ/dt|^2> 0.120310 eV Occupation number 26.756361 Fit temperature 1395.986797 K Base line 0.000553 eV * ps Maximum height 0.042598 eV * ps Fitting global error 0.025425 Frequency shift -0.255003 THz Peak # 3 ---------------------------------------------- Width 0.847729 THz Position 8.762397 THz Area (<K>) (Lorentzian) 0.009423 eV Area (<K>) (Total) 0.011871 eV <|dQ/dt|^2> 0.018846 eV Occupation number 2.767683 Fit temperature 216.985862 K Base line 0.000127 eV * ps Maximum height 0.007077 eV * ps Fitting global error 0.078128 Frequency shift -0.243450 THz Peak # 4 ---------------------------------------------- Width 0.822505 THz Position 13.296042 THz Area (<K>) (Lorentzian) 0.012441 eV Area (<K>) (Total) 0.016789 eV <|dQ/dt|^2> 0.024882 eV Occupation number 2.343156 Fit temperature 285.744350 K Base line 0.000221 eV * ps Maximum height 0.009629 eV * ps Fitting global error 0.060141 Frequency shift -0.428874 THz Peak # 5 ---------------------------------------------- Width 0.692941 THz Position 14.927237 THz Area (<K>) (Lorentzian) 0.037140 eV Area (<K>) (Total) 0.038953 eV <|dQ/dt|^2> 0.074280 eV Occupation number 7.060040 Fit temperature 860.720269 K Base line 0.000129 eV * ps Maximum height 0.034121 eV * ps Fitting global error 0.019502 Frequency shift -0.499209 THz Peak # 6 ---------------------------------------------- Width 0.873800 THz Position 15.089642 THz Area (<K>) (Lorentzian) 0.020289 eV Area (<K>) (Total) 0.022192 eV <|dQ/dt|^2> 0.040578 eV Occupation number 3.585462 Fit temperature 468.522605 K Base line 0.000120 eV * ps Maximum height 0.014782 eV * ps Fitting global error 0.026540 Frequency shift -0.493114 THz 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 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... 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... 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.531399 THz Position 4.473542 THz Area (<K>) (Lorentzian) 0.029344 eV Area (<K>) (Total) 0.027604 eV <|dQ/dt|^2> 0.058688 eV Occupation number 19.431012 Fit temperature 680.898947 K Base line -0.000048 eV * ps Maximum height 0.035154 eV * ps Fitting global error 0.015956 Frequency shift -0.194337 THz Peak # 2 ---------------------------------------------- Width 0.531399 THz Position 4.473542 THz Area (<K>) (Lorentzian) 0.029344 eV Area (<K>) (Total) 0.027604 eV <|dQ/dt|^2> 0.058688 eV Occupation number 19.431012 Fit temperature 680.898947 K Base line -0.000048 eV * ps Maximum height 0.035154 eV * ps Fitting global error 0.015956 Frequency shift -0.194337 THz Peak # 3 ---------------------------------------------- Width 0.538267 THz Position 10.993312 THz Area (<K>) (Lorentzian) 0.062518 eV Area (<K>) (Total) 0.058145 eV <|dQ/dt|^2> 0.125036 eV Occupation number 16.779904 Fit temperature 1450.579466 K Base line -0.000158 eV * ps Maximum height 0.073941 eV * ps Fitting global error 0.008088 Frequency shift -0.317902 THz Peak # 4 ---------------------------------------------- Width 0.776728 THz Position 12.855329 THz Area (<K>) (Lorentzian) 0.017087 eV Area (<K>) (Total) 0.019549 eV <|dQ/dt|^2> 0.034174 eV Occupation number 3.538695 Fit temperature 394.533162 K Base line 0.000136 eV * ps Maximum height 0.014005 eV * ps Fitting global error 0.044153 Frequency shift -0.299509 THz Peak # 5 ---------------------------------------------- Width 0.632792 THz Position 14.950742 THz Area (<K>) (Lorentzian) 0.063192 eV Area (<K>) (Total) 0.069609 eV <|dQ/dt|^2> 0.126384 eV Occupation number 12.342941 Fit temperature 1465.887345 K Base line 0.000371 eV * ps Maximum height 0.063574 eV * ps Fitting global error 0.012098 Frequency shift -0.475704 THz Peak # 6 ---------------------------------------------- Width 0.632792 THz Position 14.950742 THz Area (<K>) (Lorentzian) 0.063192 eV Area (<K>) (Total) 0.069609 eV <|dQ/dt|^2> 0.126384 eV Occupation number 12.342941 Fit temperature 1465.887345 K Base line 0.000371 eV * ps Maximum height 0.063574 eV * ps Fitting global error 0.012098 Frequency shift -0.475704 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.75 0.5 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.75 0.5 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.75 0.5 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]
Out[15]:
array([[[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+01]],
[[-8.43844610e-03, 7.27856474e-04, -8.95058148e-04],
[ 7.27856474e-04, 7.28656134e-03, 8.95058148e-04],
[-8.95058148e-04, 8.95058148e-04, 8.65518224e-03]],
[[ 7.28656134e-03, 7.27856474e-04, 8.95058148e-04],
[ 7.27856474e-04, -8.43844610e-03, -8.95058148e-04],
[ 8.95058148e-04, -8.95058148e-04, 8.65518224e-03]],
...,
[[-4.34304104e-03, 1.31864165e-04, -2.95050844e-03],
[ 1.31864165e-04, 3.82031823e-03, 9.26457597e-03],
[-2.95050844e-03, 9.26457597e-03, 2.59330404e-03]],
[[ 3.82031823e-03, 1.31864165e-04, 9.26457597e-03],
[ 1.31864165e-04, -4.34304104e-03, -2.95050844e-03],
[ 9.26457597e-03, -2.95050844e-03, 2.59330404e-03]],
[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]]],
[[[-8.43844610e-03, 7.27856474e-04, -8.95058148e-04],
[ 7.27856474e-04, 7.28656134e-03, 8.95058148e-04],
[-8.95058148e-04, 8.95058148e-04, 8.65518224e-03]],
[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+01]],
[[ 5.65686089e-03, -7.27856474e-04, -8.95058148e-04],
[-7.27856474e-04, 5.65686089e-03, -8.95058148e-04],
[-8.95058148e-04, -8.95058148e-04, 1.29289260e-02]],
...,
[[-3.48118608e+00, -2.36179399e+00, 2.36461264e+00],
[-2.36179399e+00, -3.48118608e+00, 2.36461264e+00],
[ 2.36461264e+00, 2.36461264e+00, -3.47995906e+00]],
[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]],
[[ 3.82031823e-03, 1.31864165e-04, 9.26457597e-03],
[ 1.31864165e-04, -4.34304104e-03, -2.95050844e-03],
[ 9.26457597e-03, -2.95050844e-03, 2.59330404e-03]]],
[[[ 7.28656134e-03, 7.27856474e-04, 8.95058148e-04],
[ 7.27856474e-04, -8.43844610e-03, -8.95058148e-04],
[ 8.95058148e-04, -8.95058148e-04, 8.65518224e-03]],
[[ 5.65686089e-03, -7.27856474e-04, -8.95058148e-04],
[-7.27856474e-04, 5.65686089e-03, -8.95058148e-04],
[-8.95058148e-04, -8.95058148e-04, 1.29289260e-02]],
[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+01]],
...,
[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]],
[[-3.48118608e+00, -2.36179399e+00, 2.36461264e+00],
[-2.36179399e+00, -3.48118608e+00, 2.36461264e+00],
[ 2.36461264e+00, 2.36461264e+00, -3.47995906e+00]],
[[-4.34304104e-03, 1.31864165e-04, -2.95050844e-03],
[ 1.31864165e-04, 3.82031823e-03, 9.26457597e-03],
[-2.95050844e-03, 9.26457597e-03, 2.59330404e-03]]],
...,
[[[-4.34304104e-03, 1.31864165e-04, -2.95050844e-03],
[ 1.31864165e-04, 3.82031823e-03, 9.26457597e-03],
[-2.95050844e-03, 9.26457597e-03, 2.59330404e-03]],
[[-3.48118608e+00, -2.36179399e+00, 2.36461264e+00],
[-2.36179399e+00, -3.48118608e+00, 2.36461264e+00],
[ 2.36461264e+00, 2.36461264e+00, -3.47995906e+00]],
[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]],
...,
[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+01]],
[[ 5.65686089e-03, -7.27856474e-04, -8.95058148e-04],
[-7.27856474e-04, 5.65686089e-03, -8.95058148e-04],
[-8.95058148e-04, -8.95058148e-04, 1.29289260e-02]],
[[ 7.28656134e-03, 7.27856474e-04, 8.95058148e-04],
[ 7.27856474e-04, -8.43844610e-03, -8.95058148e-04],
[ 8.95058148e-04, -8.95058148e-04, 8.65518224e-03]]],
[[[ 3.82031823e-03, 1.31864165e-04, 9.26457597e-03],
[ 1.31864165e-04, -4.34304104e-03, -2.95050844e-03],
[ 9.26457597e-03, -2.95050844e-03, 2.59330404e-03]],
[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]],
[[-3.48118608e+00, -2.36179399e+00, 2.36461264e+00],
[-2.36179399e+00, -3.48118608e+00, 2.36461264e+00],
[ 2.36461264e+00, 2.36461264e+00, -3.47995906e+00]],
...,
[[ 5.65686089e-03, -7.27856474e-04, -8.95058148e-04],
[-7.27856474e-04, 5.65686089e-03, -8.95058148e-04],
[-8.95058148e-04, -8.95058148e-04, 1.29289260e-02]],
[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+01]],
[[-8.43844610e-03, 7.27856474e-04, -8.95058148e-04],
[ 7.27856474e-04, 7.28656134e-03, 8.95058148e-04],
[-8.95058148e-04, 8.95058148e-04, 8.65518224e-03]]],
[[[ 4.40243519e-03, -7.86267323e-04, 7.46626817e-03],
[-7.86267323e-04, 4.40243519e-03, 7.46626817e-03],
[ 7.46626817e-03, 7.46626817e-03, 1.87505791e-03]],
[[ 3.82031823e-03, 1.31864165e-04, 9.26457597e-03],
[ 1.31864165e-04, -4.34304104e-03, -2.95050844e-03],
[ 9.26457597e-03, -2.95050844e-03, 2.59330404e-03]],
[[-4.34304104e-03, 1.31864165e-04, -2.95050844e-03],
[ 1.31864165e-04, 3.82031823e-03, 9.26457597e-03],
[-2.95050844e-03, 9.26457597e-03, 2.59330404e-03]],
...,
[[ 7.28656134e-03, 7.27856474e-04, 8.95058148e-04],
[ 7.27856474e-04, -8.43844610e-03, -8.95058148e-04],
[ 8.95058148e-04, -8.95058148e-04, 8.65518224e-03]],
[[-8.43844610e-03, 7.27856474e-04, -8.95058148e-04],
[ 7.27856474e-04, 7.28656134e-03, 8.95058148e-04],
[-8.95058148e-04, 8.95058148e-04, 8.65518224e-03]],
[[ 1.47901738e+01, -7.27856474e-04, 8.95058148e-04],
[-7.27856474e-04, 1.47901738e+01, 8.95058148e-04],
[ 8.95058148e-04, 8.95058148e-04, 1.47888052e+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:
In [16]:
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.
In [17]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps_library, 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'
)
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,
diable_log_file=True,
)
eng_pot_lst, eng_kin_lst = [], []
for i in range(steps):
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps_library(
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
Phonons¶
To calculate the phonons at a fixed volume the PhonopyWorkflow is used:
In [18]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, get_potential_by_name
from atomistics.workflows import PhonopyWorkflow
from phonopy.units import VaspToTHz
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
workflow = PhonopyWorkflow(
structure=bulk("Al", cubic=True),
interaction_range=10,
factor=VaspToTHz,
displacement=0.01,
dos_mesh=20,
primitive_matrix=None,
number_of_snapshots=None,
)
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
phonopy_dict = workflow.analyse_structures(output_dict=result_dict)
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 atomsfactorconversion factor, typically justphonopy.units.VaspToTHzdisplacementdisplacement 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:
In [19]:
tp_dict = workflow.get_thermal_properties(
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]), '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.10363972e-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]), '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])}
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:
In [20]:
mat = workflow.get_dynamical_matrix()
mat
Out[20]:
array(None, dtype=object)
Or alternatively the hesse matrix:
In [21]:
mat = workflow.get_hesse_matrix()
mat
Out[21]:
array([[ 4.50127147e-02, -1.92714960e-33, 8.52306995e-33, ...,
-6.63514216e-05, 8.82979633e-06, 5.93920137e-05],
[-5.07378488e-34, 4.50127147e-02, 5.07378488e-34, ...,
8.82979633e-06, -6.63514216e-05, 5.93920137e-05],
[ 5.07378488e-34, -5.07378488e-34, 4.50127147e-02, ...,
5.93659141e-05, 5.93659141e-05, 1.73512126e-05],
...,
[-6.63514216e-05, 8.82979633e-06, 5.93920137e-05, ...,
4.50127147e-02, -1.92714960e-33, 8.52306995e-33],
[ 8.82979633e-06, -6.63514216e-05, 5.93920137e-05, ...,
-5.07378488e-34, 4.50127147e-02, 5.07378488e-34],
[ 5.93659141e-05, 5.93659141e-05, 1.73512126e-05, ...,
5.07378488e-34, -5.07378488e-34, 4.50127147e-02]])
Finally, also the function to calculate the band structure is directly available on the PhonopyWorkflow:
In [22]:
band_structure = workflow.get_band_structure(
npoints=101,
with_eigenvectors=False,
with_group_velocities=False
)
This band structure can also be visualised using the built-in plotting function:
In [23]:
workflow.plot_band_structure()
Out[23]:
<Axes: title={'center': 'Bandstructure'}, xlabel='Bandpath', ylabel='Frequency [THz]'>
Just like the desnsity of states which can be plotted using:
In [24]:
workflow.plot_dos()
Out[24]:
<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:
In [25]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, get_potential_by_name
from atomistics.workflows import QuasiHarmonicWorkflow
potential_dataframe = get_potential_by_name(
potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1'
)
workflow = QuasiHarmonicWorkflow(
structure=bulk("Al", cubic=True),
num_points=11,
vol_range=0.05,
interaction_range=10,
factor=VaspToTHz,
displacement=0.01,
dos_mesh=20,
primitive_matrix=None,
number_of_snapshots=None,
)
task_dict = workflow.generate_structures()
result_dict = evaluate_with_lammps(
task_dict=task_dict,
potential_dataframe=potential_dataframe,
)
fit_dict = workflow.analyse_structures(output_dict=result_dict)
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:
In [26]:
tp_dict = workflow.get_thermal_properties(
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]), 'free_energy': array([ 0.1490366 , 0.14826761, 0.13934538, 0.11699115, 0.08193384,
0.03597288, -0.01929866, -0.08261783, -0.15299316, -0.22963661,
-0.31191108, -0.39929276, -0.49134426, -0.58769537, -0.68802891,
-0.7920703 , -0.89957955, -1.01034525, -1.12417973, -1.24091537,
-1.36040152, -1.48250213, -1.6070937 , -1.73406366, -1.86330899,
-1.99473507, -2.12825469, -2.26378725, -2.40125802, -2.54059753,
-2.68174106]), 'entropy': array([1.36944221e-05, 5.98139595e+00, 2.96871756e+01, 5.55859383e+01,
7.82416454e+01, 9.75225732e+01, 1.14066300e+02, 1.28465948e+02,
1.41175404e+02, 1.52531113e+02, 1.62783732e+02, 1.72122876e+02,
1.80694518e+02, 1.88612990e+02, 1.95969277e+02, 2.02836853e+02,
2.09275827e+02, 2.15335964e+02, 2.21058900e+02, 2.26479810e+02,
2.31628669e+02, 2.36531221e+02, 2.41209738e+02, 2.45683613e+02,
2.49969835e+02, 2.54083369e+02, 2.58037467e+02, 2.61843913e+02,
2.65513240e+02, 2.69054894e+02, 2.72477381e+02]), 'heat_capacity': array([6.93153576e-05, 1.73540236e+01, 5.38037703e+01, 7.36871467e+01,
8.35644373e+01, 8.88841670e+01, 9.20085314e+01, 9.39792226e+01,
9.52946943e+01, 9.62133949e+01, 9.68788949e+01, 9.73756860e+01,
9.77559502e+01, 9.80532531e+01, 9.82899461e+01, 9.84813617e+01,
9.86382928e+01, 9.87685099e+01, 9.88777194e+01, 9.89701870e+01,
9.90491514e+01, 9.91171071e+01, 9.91759998e+01, 9.92273650e+01,
9.92724271e+01, 9.93121723e+01, 9.93474015e+01, 9.93787710e+01,
9.94068222e+01, 9.94320052e+01, 9.94546965e+01]), 'volumes': [66.7171076313667, 66.72177243343418, 66.75880349778612, 66.82252096053044, 66.89849531395458, 66.9796341067605, 67.06260761039168, 67.14571326277367, 67.22800273919441, 67.30891229729374, 67.38809258770314, 67.46532342705305, 67.5404676964542, 67.61344459222596, 67.6842130824627, 67.75276107005782, 67.81909792596628, 67.88324911963028, 67.94525222117329, 68.00515384392953, 68.06300725960926, 68.11887051273405, 68.17280491720088, 68.22487385247723, 68.27514179905809, 68.32367356749006, 68.37053368537661, 68.41578591403304, 68.45949287183225, 68.50171574542424, 68.54251407326704]}
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:
In [27]:
temperatures, volumes = workflow.get_thermal_expansion(
output_dict=result_dict,
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,
)
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:
In [28]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, 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'
)
result_dict = evaluate_with_lammps(
task_dict=optimize_positions_and_volume(structure=structure),
potential_dataframe=potential_dataframe,
)
structure_opt = result_dict["structure_with_optimized_positions_and_volume"]
structure_opt
Out[28]:
Atoms(symbols='Al4', pbc=True, cell=[[4.05000466219724, 2.4799126230458533e-16, 2.4799126230458533e-16], [0.0, 4.05000466219724, 2.4799126230458533e-16], [0.0, 0.0, 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:
In [29]:
from ase.build import bulk
from atomistics.calculators import evaluate_with_lammps, 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'
)
result_dict = evaluate_with_lammps(
task_dict=optimize_positions(structure=structure),
potential_dataframe=potential_dataframe,
)
structure_opt = result_dict["structure_with_optimized_positions"]
structure_opt
Out[29]:
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.
In [ ]: