DynamicalMatrix¶

class
DynamicalMatrix
(configuration, filename, object_id, calculator=None, repetitions=None, atomic_displacement=None, acoustic_sum_rule=None, finite_difference_method=None, constraints=None, constrain_electrodes=None, use_equivalent_bulk=None, max_interaction_range=None, force_tolerance=None, processes_per_displacement=1, log_filename_prefix='forces_displacement_', use_wigner_seitz_scheme=None, use_symmetry=None, symmetrize=None)¶ Constructor for the DynamicalMatrix object.
Parameters:  configuration (
BulkConfiguration
MoleculeConfiguration
DeviceConfiguration
) – The configuration for which to calculate the dynamical matrix.  calculator (
Calculators
) – The calculator to be used in the dynamical matrix calculations.
Default: The calculator attached to the configuration.  filename (str) – The full or relative path to save the results to. See
nlsave()
.  object_id (str) – The object id to use when saving. See
nlsave()
.  repetitions (
Automatic
 list of ints) – The number of repetitions of the system in the A, B, and Cdirections given as a list of three positive integers, e.g.[3, 3, 3]
, orAutomatic
. Each repetition value must be odd. Ifuse_wigner_seitz_scheme
is set toTrue
the only values allowed for repetition are[1, 1, 1]
orAutomatic
.
Default:Automatic
 atomic_displacement (PhysicalQuantity of type length) – The distance the atoms are displaced in the finite difference method.
Default:0.01 * Angstrom
 acoustic_sum_rule (bool) – Control if the acoustic sum rule should be invoked.
Default:True
 finite_difference_method (
Forward
Central
) – The finite difference scheme to use.
Default:Central
 constraints (list of type int) – List of atomic indices that will be constrained, e.g.
[0, 2, 10]
.
Default: Empty list[]
 constrain_electrodes (bool) – Control if the electrodes and electrode extensions should be constrained in case of a
DeviceConfiguration
.
Default:False
 use_equivalent_bulk (bool) – Control if a
DeviceConfiguration
should be treated as aBulkConfiguration
.
Default:True
 max_interaction_range (PhysicalQuantity of type length) – Set the maximum range of the interactions.
Default: All atoms are included  force_tolerance (PhysicalQuantity of type energy per length squared) – All force constants below this value will be truncated to zero.
Default:1e8 * Hartree/Bohr**2
 processes_per_displacement (int) – The number of processes assigned to calculating a single displacement.
Default: 1 process per displacement.  log_filename_prefix (str or None) – Prefix for the filenames where the logging output for every displacement calculation is
stored. The filenames are formed by appending a number and the file extension (”.log”).
If a value of None is given then all logging output is done to stdout. If a classical
calculator is used, no perdisplacment log files will be generated.
Default:"displacement_"
 use_wigner_seitz_scheme (bool) – Control if the real space Dynamical Matrix should be extended according to the Wigner
Seitz construction.
use_wigner_seitz_scheme=True
is only supported for simple orthorhombic, simple tetragonal and simple cubic lattices.
Default:False
 use_symmetry (bool) – If enabled, only the symmetrically unique atoms are displaced and the remaining force
constants are calculated using symmetry.
Default:True

acousticSumRule
()¶ Returns: Return if the acoustic sum rule is invoked. Return type: bool

atomicDisplacement
()¶ Returns: The distance the atoms are displaced in the finite difference method. Return type: PhysicalQuantity with length unit

constrainElectrodes
()¶ Returns: Boolean determining if the electrodes and electrode extensions are constrained in case of a DeviceConfiguration. Return type: bool

constraints
()¶ Returns: The list of constrained atoms. Return type: list of int

filename
()¶ Returns: The filename where the study object is stored. Return type: str

finiteDifferenceMethod
()¶ Returns: The finite difference scheme to use. Return type: Central
Forward

forceTolerance
()¶ Returns: The force tolerance Return type: PhysicalQuantity with an energy per length squared unit e.g. Hartree/Bohr**2

logFilenamePrefix
()¶ Returns: The filename prefix for the logging output of the study. Return type: str  LogToStdOut

maxInteractionRange
()¶ Returns: The maximum interaction range. Return type: PhysicalQuantity with length unit

nlprint
(stream=None)¶ Print a string containing an ASCII table useful for plotting the Study object.
Parameters: stream (python stream) – The stream the table should be written to.
Default:NLPrintLogger()

numberOfProcessesPerTask
()¶ Returns: The number of processes to be used to execute each task. If None, all available processes should execute each task collaboratively. Return type: int  None

objectId
()¶ Returns: The name of the study object in the file. Return type: str

phononEigensystem
(q_point=None, constrained_atoms=None)¶ Calculate the eigenvalues and eigenvectors for the dynamical matrix at a specified qpoint.
Parameters:  q_point (list of 3 floats) – The fractional qpoint to use.
Default:[0.0, 0.0, 0.0]
 constrained_atoms (list of ints.) – List of atoms being constrained. The matrix elements from these atoms will not be
included in the calculation of the eigensystem.
Default: [] (empty list)
Returns: The eigenvalues and eigenvectors of the dynamical matrix.
Return type: 2tuple containing the eigenvalues and eigenvectors
 q_point (list of 3 floats) – The fractional qpoint to use.

processesPerDisplacement
()¶ Returns: The number of processes per displacement. Return type: int

realSpaceDynamicalMatrix
()¶ Returns the real space dynamical matrix. The shape of the matrix is (N, M), where N is the number of degrees of freedom (3 * number of atoms), and M = N * R, where R is the total number of repetitions.
Each subblock D[i*N:(i+1)*N, :] corresponds to the matrix elements between the center block (where the atoms have been displaced) and a neighbouring cell translated from the central cell by translations[i] in fractional coordinates.
Returns: The realspace dynamical matrix as a sparse matrix together with a list of translation vectors in fractional coordinates. The real space dynamical matrix is given in units of (meV / hbar)**2. Return type: (scipy.sparse.csr_matrix, list of list(3) of integers)

reciprocalSpaceDynamicalMatrix
(q_point=None)¶ Evaluate the reciprocal space dynamical matrix for a given qpoint in reciprocal space.
Parameters: q_point (list of floats) – The fractional qpoint to use.
Default:[0.0, 0.0, 0.0]
Returns: The dynamical matrix for q_point
.Return type: PhysicalQuantity of units (meV / hbar)**2

repetitions
()¶ Returns: The number of repetitions in the A, B, and Cdirections for the supercell that is used in the finite displacement calculation. Return type: list of three int.

symmetry
()¶ Returns: True if the use of crystal symmetry to reduce the number of displacements is enabled. Return type: bool

update
()¶ Run the calculations for the DynamicalMatrix study object.

useEquivalentBulk
()¶ Returns: Boolean determining if a DeviceConfiguration is treated as a BulkConfiguration. Return type: bool

wignerSeitzScheme
()¶ Returns: Boolean to control if the real space Dynamical Matrix should be extended according to the WignerSeitz construction. Return type: bool
 configuration (
Usage Examples¶
Note
Study objects behave differently from analysis objects. See the Study object overview for more details.
Calculate the DynamicalMatrix
for a system repeated five times in the B direction and three times in the C direction.
dynamical_matrix = DynamicalMatrix(
configuration,
filename='DynamicalMatrix.hdf5',
object_id='dynamical_matrix',
repetitions=(1,5,3)
)
dynamical_matrix.update()
When using repetitions=Automatic
, the cell is repeated such that all atoms within a predefined, elementpair dependent interaction range are included.
dynamical_matrix = DynamicalMatrix(
configuration,
filename='DynamicalMatrix.hdf5',
object_id='dynamical_matrix',
repetitions=Automatic
)
dynamical_matrix.update()
The default number of repetitions i.e. repetitions=Automatic
can be found before a calculation using the function
checkNumberOfRepetitions()
.
(nA, nB, nC) = checkNumberOfRepetitions(configuration)
The maximum interaction range between two atoms can be specified manually, by using
the max_interaction_range
keyword.
dynamical_matrix = DynamicalMatrix(
configuration,
filename='DynamicalMatrix.hdf5',
object_id='dynamical_matrix',
repetitions=Automatic,
max_interaction_range=12.0*Ang,
)
dynamical_matrix.update()
Notes¶
The DynamicalMatrix
is calculated using the finite difference
method in a repeated cell, which is sometimes also referred to as frozen
phonon or supercell method.
In the following, we denote the atoms in the unit cell by \(\mu\) and the atoms in the repeated cell by \(i\). Furthermore, denote the dynamical matrix elements, \(D_{\mu \alpha, i \beta}\), where \(\alpha, \beta\) are the Cartesian directions, i.e. \(x, y, z\). A dynamical matrix element is given by
where \(F_{i \beta}\) is the force on atom \(i\) in direction \(\beta\) due to a displacement of atom \(\mu\) in direction \(\alpha\).
The derivative is calculated by either forward
or central
finite
differences, where in the following we will focus on the latter. Atom
\(\mu\) is displaced by \(\Delta r_\alpha\) and \(\Delta
r_\alpha\), and the changes in the force, \(\Delta F_{i \beta}\) are
calculated to approximate the dynamical matrix element
The default is to use repetitions=Automatic
. In this case the cell is
repeated such that all atoms within a predefined elementdependent distance
from the atoms in the unit cell are included in the repeated cell. The
repetitions used is written to the output file. These default interaction
ranges are suitable for most systems.
However, if you are using longranged interactions, e.g. classical potentials
with electrostatic interactions in TremoloXCalculator
, it might be
necessary to increase the number of repetitions.
For a 1D or 2D system, the unit cell should not be repeated in the confined
directions. This is only discovered by the repetitions=Automatic
, if there
is enough vacuum in the unit cell in the directions that should not be
repeated. That is typically 1020 Å vacuum depending on the elements and their
interaction ranges. Thus, for confined systems it is recommended to check the
repetitions used and possibly use manual instead of automatic repetition.
DynamicalMatrix calculations using DFT or SemiEmpirical calculators have
functionality to fully resume partially completed calculations by rerunning the
same script or reading the study object from file and calling update()
on
it. The study object will automatically detect which displacement
calculations have already been carried out and only run the ones that are not
yet completed. To ensure highest performance this functionality is not available
for ATKForceField calculations.
Notes for DFT¶
In ATKDFT the number of repetitions of the unit cell in super cell must ensure
that the change in the force on atoms outside the super cell is zero for every
atomic displacement in the center cell. An equivalent discussion of the number
of kpoints of the super cell and the number of repetitions can be found for
HamiltonianDerivatives
in the section Notes for DFT.
Consider a system with e.g. \(x\) and \(y\) as confined directions and
the kpoint sampling of the unit cell \((1, 1, N_{\text{C}})\), see
Fig. 131 (a). Assume that the number of
repetitions in the Cdirection is known for the change in the force on atoms
outside the super cell to be zero. Then the number of repetitions must be
\((1, 1, {\scriptstyle \text{repetitions in C}})\). Furthermore, the kpoint
sampling of the super cell becomes \((1, 1, \frac{N_{\text{C}}}{\text{repetitions in C}})\).
Note
From QuantumATK2019.03 onwards, the kpoint sampling and
densitymeshcutoff will be automatically adapted to the given number
of repetitions when setting up the super cell inside
DynamicalMatrix
and HamiltonianDerivatives
. That
means you can specify the calculator settings for the unit cell and use it
with any desired number of repetitions in dynamical matrix and hamiltonian
derivatives calculations.
When calculating the DynamicalMatrix
with ATKDFT, accurate results
may require a higher precision than usual by increasing the density_mesh_cutoff
in NumericalAccuracyParameters
and decreasing the tolerance
in
IterationControlParameters
, e.g.
numerical_accuracy_parameters = NumericalAccuracyParameters(
density_mesh_cutoff=150.0*Hartree
)
iteration_control_parameters = IterationControlParameters(
tolerance=1e6
)
Notes for the simplified supercell method (use_wigner_seitz_scheme=True
)¶
The simplified supercell method is an approximation which allows to obtain the
dispersion of vibrational eigenmodes with a force calculation of the unit cell
only, i.e. having repetitions=[1,1,1]
(the poor man’s frozen phonon
calculation). It is valid if the unit cell is large enough, i.e. if it
accommodates 200 atoms and more. The convergence should be checked with
respect to the number of atoms per unit cell or by a conventional frozen
phonon calculation.
Idea of the simplified supercell method¶
In large unit cells, the force of atom i due to a displacement of atom j is small for distances of half the length of the unit cell vectors. Due to the translational invariance, a displaced atom i contributes from two sides to the force on atom j, which can be exactly decomposed if the distance between the atoms is large enough. As an example, we will look at a 1D chain with 6 atoms per unit cell (see Fig. 127).
All atoms are at their equilibrium positions, besides atom 1 which is
displaced along the \(x\) direction. The force on atom 5, for example,
can be regarded as “repulsive” due to the displacement of atom 1 in the unit
cell with translation \(T=0\) along the \(x\)axis. However, there is
also an “attractive” contribution from the image of the displaced atom 1 in
the unit cell with translation \(T=1\). To decompose the two contributions
we construct a WignerSeitz cell centered at the undisplaced position of atom
1 and check the distance to the nearest neighbor representations of atom 5.
The representation of atom 5 inside the cell at \(T=0\) is not part of the
WignerSeitz cell, and thus this contribution is neglected. In contrast, the
representation of atom 5 at \(T=1\) is inside the WignerSeitz cell and
we keep this contribution. If an atom is at the face (like atom 4), the edge
or a vertex of the WignerSeitz cell, the force is included weighted by the
inverse multiplicity of this WignerSeitz grid point. This construction is
repeated for all atoms inside the unit cell. In this way it is possible to
approximately calculate the entries of the DynamicalMatrix
for all
repetitions of the cell from \(T = 1\) to \(T = 1\) and thereby get
the dispersion, by only calculating the forces in a single unit cell.
Note
The WignerSeitz construction ensures that the phonon spectrum at the \(\Gamma\)point is preserved. Hence, a separate \(\Gamma\)point calculation is not necessary.