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 C-directions given as a list of three positive integers, e.g. [3, 3, 3], or Automatic. Each repetition value must be odd. If use_wigner_seitz_scheme is set to True the only values allowed for repetition are [1, 1, 1] or Automatic.
    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 a BulkConfiguration.
    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: 1e-8 * 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 per-displacment 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 q-point.

Parameters:
  • q_point (list of 3 floats) – The fractional q-point 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:

2-tuple containing the eigenvalues and eigenvectors
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 real-space 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 q-point in reciprocal space.

Parameters:q_point (list of floats) – The fractional q-point 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 C-directions 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 Wigner-Seitz construction.
Return type:bool

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 pre-defined, element-pair 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 super-cell 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

\[D_{\mu \alpha, i \beta} = \frac{d F_{i \beta}}{d r_{\mu \alpha}},\]

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

\[D_{\mu \alpha, i \beta} \approx \frac{F_{i \beta}(\Delta r_\alpha) - F_{i \beta}(-\Delta r_\alpha)}{2\Delta r_\alpha}.\]

The default is to use repetitions=Automatic. In this case the cell is repeated such that all atoms within a pre-defined element-dependent 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 long-ranged 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 10-20 Å 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 Semi-Empirical calculators have functionality to fully resume partially completed calculations by re-running 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 ATK-ForceField calculations.

Notes for DFT

In ATK-DFT 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 k-points 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 k-point sampling of the unit cell \((1, 1, N_{\text{C}})\), see Fig. 131 (a). Assume that the number of repetitions in the C-direction 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 k-point sampling of the super cell becomes \((1, 1, \frac{N_{\text{C}}}{\text{repetitions in C}})\).

Note

From QuantumATK-2019.03 onwards, the k-point sampling and density-mesh-cutoff 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 ATK-DFT, 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=1e-6
    )

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).

../../../_images/simplified_supercell_method_explained.png

Fig. 127 A 1D chain with 6 atoms per unit cell with atom 1 being displaced. The periodically-repeated atoms to the left (\(T=-1\)) and to the right (\(T=1\)) are indicated. The Wigner-Seitz cell centered at the undisplaced position of atom 1 is drawn red / dashed.

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 Wigner-Seitz 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 Wigner-Seitz cell, and thus this contribution is neglected. In contrast, the representation of atom 5 at \(T=-1\) is inside the Wigner-Seitz cell and we keep this contribution. If an atom is at the face (like atom 4), the edge or a vertex of the Wigner-Seitz cell, the force is included weighted by the inverse multiplicity of this Wigner-Seitz 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 Wigner-Seitz construction ensures that the phonon spectrum at the \(\Gamma\)-point is preserved. Hence, a separate \(\Gamma\)-point calculation is not necessary.