SOFI reference

Functions

subroutine sofi_compute_all(nat, typ, coords, sym_thr, prescreen_ih, nmat, mat_list, perm_list, op_list, n_list, p_list, ax_list, angle_list, dHausdorff_list, pg, n_prin_ax, prin_ax, ierr)

Compute all data from SOFI in single routine. This contains getting the list of symmetry operations, and the associated permutations, Op symbols, n and p integers, list of axis, angles, and dHausdorff. Also the PG name, and list of principal axes.

The structure is assumed to be already shifted to desired origin on input.

Note

All arrays have size nmax, but the actual output values go only up to index nmat (or n_prin_ax), values beyond that index can be random.

Warning

This routine performs the combinations of found symmetry operations until group completeness. The sym_thr is disregarded for the combinations, therefore operations resulting as combinations of other operations can have a dHausdorff value higher than sym_thr!

Parameters

nat – [in] :: number of atoms
typ(nat) – [in] :: integer atomic types
coords(3, nat) – [in] :: positions of atoms
sym_thr – [in] :: threshold for finding symmetries (not taken into account when making combinations)
prescreen_ih – [in] :: flag to check early-termination for Ih groups
nmat – [out] :: number of symmetries found
mat_list(3, 3, nmat) – [out] :: symmetry matrices
perm_list(nat, nmat) – [out] :: permutation of atoms after applying each symmetry
op_list(nmat) – [out] :: Character “Op” from the Schoenflies notation: Op n^p (E = identity, I = inversion, C = rotation, S = (roto-)reflection )
n_list(nmat) – [out] :: Schoenflies n value
p_list(nmat) – [out] :: Schoenflies p value
ax_list(3, nmat) – [out] :: axis of operation of each symmetry
angle_list(nmat) – [out] :: angle of each symmetry, in units of 1/2pi, i.e. angle=0.333 is 1/3 of full circle, or 120 degrees
dHausdorff_list(nmat) – [out] :: max difference of atomic positions of before/after symm transformation
pg – [out] :: name of Point group, e.g. D6h
n_prin_ax – [out] :: number of equivalent principal axes
prin_ax(3, n_prin_ax) – [out] :: list of equivalent principal axes.
ierr – [out] :: error value, zero on normal execution, negative otherwise

Returns

nmat, mat_list, perm_list, op_list, n_list, p_list, ax_list, angle_list, dHausdorff_list, pg, prin_ax, ierr

subroutine sofi_struc_pg(nat, typ_in, coords_in, sym_thr, pg, verb)

Wrapper routine for directly getting only the PG name from a structure, nothing else. The structure is shifted to geometrical center in this routine!

Warning

This routine performs the combinations of found symmetry operations until group completeness. The sym_thr is disregarded for the combinations, therefore operations resulting as combinations of other operations can have a dHausdorff value higher than sym_thr!

Parameters

nat – [in] :: number of atoms
typ_in(nat) – [in] :: atomic species
coords_in(3, nat) – [in] :: atomic positions
sym_thr – [in] :: symmetry threshold
pg – [out] :: point group
verb – [in] :: flag for verbose output

subroutine sofi_get_symmops(nat, typ_in, coords_in, sym_thr, prescreen_ih, n_so, op_list, ierr)

The main SOFI routine

Find the list of symmetry operations of the atomic structure, with a threshold sym_thr. The symmetry operations are in form of 3x3 matrices.

When distortions are present in the atomic structure, the distance between the original structure \(A\) and the symmetry-transformed structure \(\theta A\) is some small nonzero value. This distance is computed as largest distance of any atom after the transformation from the position of the corresponding atom in the original structure (Hausdorff distance). The threshold sym_thr gives the upper limit for this distance value.

The operations from op_list can be applied by matmul(), for example the N-th entry:

theta(3,3) = op_list(:, :, n)
do i = 1, natoms
   coords_transformed(:,i) = matmul( theta, coords_original(:,i) )
end do

Note

The output op_list is assumed to be allocated by the caller, and has the size [3,3,nmax], where nmax=200 by default, but can be adjusted in sofi_tools.f90

Parameters

nat – [in] :: number of atoms
typ_in(nat) – [in] :: integer atomic types
coords_in(3, nat) – [in] :: positions of atoms
sym_thr – [in] :: threshold for finding symmetries, in terms of Hausdorff distance
prescreen_ih – [in] :: flag to check early termniation for Ih
n_so – [out] :: number of found symmetry operations
op_list – [out] :: the list of operations
ierr – [out] :: error value, negative on error, zero otherwise

Returns

n_so, op_list, ierr

subroutine sofi_get_perm(nat, typ, coords, nbas, bas_list, perm_list, dHausdorff_list)

Obtain the permutations of atoms associated to each of the symmetry operations in the list. Also compute the dHausdorff of each operation.

This consists of computing P by cshda, permuting, applying SVD, and finally computing dHausdorff.

Applying the permutation associated to N-th symmetry operation:

coords(:, perm_list(:, n) ) = coords(:,:)

Parameters

nat – [in] :: number of atoms
typ(nat) – [in] :: atomic types
coords(3, nat) – [in] :: atomic positions
nbas – [in] :: number of operations in the list
bas_list(3, 3, nbas) – [in] :: list of symmetry operations
perm_list(nat, nbas) – [out] :: list of permutations
dHausdorff_list(nbas) – [out] :: list of dHausdorff

Returns

perm_list, dHausdorff_list

subroutine sofi_get_combos(nat, typ, coords, nbas, bas_list, ierr): find all combos with high sym_thr: this will accept anything. However, some strange matrices cannot be constructed from combos of matrices which are already found as symmetries with sym_thr value.

subroutine try_sofi(theta, nat, typ_in, coords_in, sym_thr, dd, nbas, op_list, dh, m_thr, success, ierr)

Try if the matrix theta gives dH within 5*sym_thr, then send it to refine. If matrix is valid and new, add it to op_list

nbas and op_list on input contain the currently known symmetry operations, and on output the same info updated accordingly if theta is new or not.

Parameters

theta – [inout] :: 3x3 trial matrix
nat – [in] :: number of atoms
typ_in(nat) – [in] :: atomic species
coords_in(3, nat) – [in] :: atomic positions
sym_thr – [in] :: symmery threshold in terms of dH
dd – [in] :: thr for first cshda, “smallest atom-atom distance”; also used for detemining “good-enough” trial matrix theta
nbas – [inout] :: number of operations
op_list(3, 3, nbas) – [inout] :: list of operations
dh – [out] :: Hausdorff distance value
m_thr – [in] :: threshold for matrix_distance checking if matrix is new or not
ierr – [out] :: error value, negative on error, zero otherwise

Returns

theta, nbas, op_list, dh

subroutine refine_sofi(nat_ref, typ_ref, coords_ref, nat, typ_in, coords_in, theta, dh, perm)

Refine the theta matrix, which is on input “close-to” some real symmetry operation, and on output should be the matrix of symmetry operation that minimizes SVD.

Do few steps of ICP-like algo… until n steps, or until no new permutations are found.

*_ref is static, original structure, *_in is the structure transformed by theta and P on input

Note

assume atoms are already permuted at input

Parameters

nat_ref – [in] :: number of atoms
typ_ref(nat_ref) – [in] :: atomic types
coords_ref(3, nat_ref) – [in] :: atomic positions
nat – [in] :: number of atoms
typ_in(nat) – [in] :: atomic types
coords_in(3, nat) – [in] :: atomic positions
theta(3, 3) – [inout] :: trial matrix
dh – [inout] :: Hausdorff value
perm(nat) – [out] :: final permutation

Returns

theta, dh, perm

subroutine is_new_sofi(rmat, nbas, op_list, m_thr, is_new)

Check if rmat is new in op_list. The choice is made by computing matrix_distance between rmat and all matrices in op_list, if the distance is below m_thr for any matrix, then rmat is considered already known.

Parameters

rmat(3, 3) – [in] :: trial matrix
nbas – [in] :: number of operations in the list
op_list(3, 3, nmax) – [in] :: list of symmetry operations to check
m_thr – [in] :: threshold for matrix_distance
is_new – [out] :: true if rmat is new, and should be added to op_list

Returns

is_new

subroutine add_sofi(nat, rmat, nbas, mat_list, ierr): Add rmat into mat_list, increase nbas by 1.

subroutine sofi_get_pg(nbas, op_list, pg, n_prin_ax, prin_ax, verb, ierr)

flowchart to determine PG notation from op_list online: https://symotter.org/assets/flowchart.pdf

Note

The principal axis is found as axis of the proper rotation operation (C symbol), with the largest n value. In case of multiple equivalent principal axes, they are found as all unique axes with those properties.

Note

In cases when the list of operations is incomplete for a certain PG, this routine will append either + or - sign to the PG tag, signifying either (+) a group with additional operations, or (-) a PG with some missing operations. These are deduced from the expected total number of operations of a PG. Such case typically happens when atomic structure is highly disordered.

Parameters

nbas – [in] :: number of symmetry operations
op_list(3, 3, nbas) – [in] :: the list of symmetry operations (3x3 matrices)
pg(10) – [out] :: the point group tag
n_prin_ax – [out] :: number of equivalent principal axes
prin_ax(3, n_prin_ax) – [out] :: list of equivalent principal axes of PG
verb – [in] :: flag for verbosity
ierr – [out] :: error value, zero on normal execution, negative otherise

Returns

pg, prin_ax

subroutine sofi_analmat(rmat, op, n, p, ax, angle, ierr)

Analyse the input 3x3 matrix rmat, return the Schoenflies PG notation: “Op n^p”, also give axis, angle of the operation.

The matrices \(M\) corresponding to PG symmetry operations are orthonormal ( \(M^{-1}=M^T\)), non-symmetric matrices with real elements, and correspond to either pure rotations (symbol \(C\)), pure reflections (symbol \(\sigma\)), or combinations of both (rotoreflections, symbol \(S\)).

The matrix properties are:

for C matrices: \(det=1\), \(\lambda_1 = 1\), and \(\lambda_{2,3}=\exp{(\pm i \alpha)}\),
for \(\sigma\) matrices: \(det=-1\), \(\lambda_1=-1\), and \(\lambda_{2,3}=1\).
for S matrices: \(det=-1\), \(\lambda_1=-1\), and \(\lambda_{2,3}=\exp{(\pm i\alpha)}\).

The axis always corresponds to the eigenvector of \(\lambda_1\).

The approach to characterise the matrices used in this routine is based on value of their determinant, eigenvalues, and eigenvectors.

The angle \(\alpha\) in radians is transformed to integers \(n\) and \(p\), such that \({\alpha}/{2\pi} = {p}/{n}\).

The use of \(\sigma\) symbol is redundant, since an equivalent notation is a rotoreflection \(S\) with angle \(\alpha=0\). We thus label pure reflections as \(S~0\).

The identity E and point-inversion I matrices always have the same form, E=diag(1,1,1) and I=diag(-1,-1,-1). Also, E is equivalent to C with angle zero, and I is equivalent to S with angle 0.5, regardless of the axis.

Note

The convention for axis oriontation is established due to the reason that two matrices M and M^T can produce the same result from the diagonalisation procedure in this routine, meaning the relative orioentation (positive/negative angle, axis direction) can be ambiguous.

Note

A different method (potentially faster) to get the angles could be based on the trace of matrix, which should be 3 for E, -3 for I, 1 for s, 2cos(a)+1 for C, 2cos(a)-1 for S. The axis can be obtained by summing all transformations of a generic vector with all matrices of the same operation and order (all p from 1 to n, or 2n for S_odd). Multiply by det(M) for S cases.

Warning

Axis output from this routine follows the convention: (1) flip ax such that z>0 (2) if z==0, then flip such that x>0 (3) if x==0, then flip such that y>0. The value of angle is then relative to this axis convention.

Parameters

rmat(3, 3) – [in] :: the input matrix
op(1) – [out] :: character for operation E, I, C, S
n – [out] :: the order of operation
p – [out] :: the power for C5^2 kinda things
ax(3) – [out] :: the axis of operation, according to convention
angle – [out] :: angle in units of 1/2pi, e.g. value 0.5 is half the circle
ierr – [out] :: error value, zero on normal execution, negative otherwise

Returns

op, n, p, ax, angle, ierr

subroutine sofi_ext_bfield(n_op, op_list, b_field)

Impose the action of an external magnetic field to the symmetry elements of the PG.

The effect of external B field on PG is to filter the list of Ops, only those which satisfy:

det( M ) M B = B

are valid, where M is the symmOp matrix.

See Eq(2) of: A. Pausch, M Gebele, W. Klopper, J. Chem. Phys. 155, 201101 (2021) https://doi.org/10.1063/5.0069859

Parameters

n_op – [inout] :: number of symmetry operations
op_list(3, 3, n_op) – [inout] :: the list of symmetry operations
b_field(3) – [in] :: direction of the B field

Returns

n_op, op_list

subroutine sofi_construct_operation(op, axis, angle, matrix, ierr)

construct a 3x3 matrix from input data of Op, axis, and angle.

A pure rotation matrix \(C\) can be constructed from knowing the angle and axis of rotation, for example by the well-known Euler-Rodrigues rotation formula. See: https://www.wikipedia/rotation_matrix#quaternion

A reflection matrix \(\sigma\) can be constructed from knowing the normal vector of the plane of reflection \(\ket{x}\) by \(\sigma=E - 2\ket{x}\bra{x}\), where \(E\) is the identity matrix, and \(\ket{.}\bra{.}\) indicates the outer product.

Matrices corresponding to rotoreflections \(S\) can be constructed from combination of pure rotation and pure reflection matrices, \(S = \sigma C\), where \(\sigma\) and \(C\) have the same vector as the axis of rotation, and the normal of the plane of reflection.

Parameters

op(1) – [in] :: character of symmetry operation E, I, C, S
axis(3) – [in] :: the axis
angle – [in] :: the angle in units 1/(2pi), e.g. angle=0.5 means half circle
matrix(3, 3) – [out] :: the orthonormal matrix
ierr – [out] :: error value, negative on error, zero otherwise

subroutine sofi_unique_ax_angle(n_mat, mat_list, op_out, ax_out, angle_out, ierr)

subroutine sofi_mat_combos(n_in, mat_in, n_out, mat_out)

subroutine sofi_get_err_msg(ierr, msg)

extract the error message from ira/sofi corresponding to error value ierr

Parameters

ierr – [in] :: integer error value
msg(128) – [out] :: error message

subroutine sofi_check_collinear(nat, coords, collinear, ax_o)

check if the atomic positions in coords form a linear structure.

Parameters

integer(ip) – [in] :: nat –> number of atoms
real(3, nat) – [in] :: coords –> atomic positions
logical – [out] :: collinear –> .true. when structure is collinear
real(3) – [out] :: ax –> axis of collinearity