Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: MIT */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : /*
9 : * libgrpp - a library for the evaluation of integrals over
10 : * generalized relativistic pseudopotentials.
11 : *
12 : * Copyright (C) 2021-2023 Alexander Oleynichenko
13 : */
14 :
15 : /*
16 : * Functions for construction of matrices of the angular momentum operator L
17 : * in the bases of either real or complex spherical harmonics.
18 : */
19 :
20 : #include "grpp_lmatrix.h"
21 :
22 : #include <complex.h>
23 : #include <math.h>
24 : #include <stdlib.h>
25 : #include <string.h>
26 :
27 : static void get_transformation_coeffs_csh_to_rsh(int m, double complex *a,
28 : double complex *b);
29 :
30 : /**
31 : * Constructs matrices of the Lx, Ly, Lz operators for the given angular
32 : * momentum L in the basis of real spherical harmonics (rsh).
33 : */
34 0 : void libgrpp_construct_angular_momentum_matrices_rsh(int L, double *lx_matrix,
35 : double *ly_matrix,
36 : double *lz_matrix) {
37 0 : int dim = 2 * L + 1;
38 :
39 : // set all matrices to zero
40 0 : memset(lx_matrix, 0, dim * dim * sizeof(double));
41 0 : memset(ly_matrix, 0, dim * dim * sizeof(double));
42 0 : memset(lz_matrix, 0, dim * dim * sizeof(double));
43 :
44 0 : double *lx_matrix_csh = calloc(dim * dim, sizeof(double));
45 0 : double *ly_matrix_csh = calloc(dim * dim, sizeof(double));
46 0 : double *lz_matrix_csh = calloc(dim * dim, sizeof(double));
47 :
48 0 : libgrpp_construct_angular_momentum_matrices_csh(L, lx_matrix_csh,
49 : ly_matrix_csh, lz_matrix_csh);
50 :
51 0 : for (int m1 = -L; m1 <= L; m1++) {
52 0 : for (int m2 = -L; m2 <= L; m2++) {
53 :
54 : // coefficients: S_lm = a * Y_{l,-m} + b * Y_{l,m}
55 0 : double complex a1, b1;
56 0 : double complex a2, b2;
57 0 : get_transformation_coeffs_csh_to_rsh(m1, &a1, &b1); // bra
58 0 : get_transformation_coeffs_csh_to_rsh(m2, &a2, &b2); // ket
59 :
60 0 : int m1m = -abs(m1);
61 0 : int m1p = +abs(m1);
62 0 : int m2m = -abs(m2);
63 0 : int m2p = +abs(m2);
64 :
65 0 : double complex lx = 0.0 + 0.0 * I;
66 0 : lx += conj(a1) * a2 * lx_matrix_csh[dim * (m1m + L) + (m2m + L)];
67 0 : lx += conj(a1) * b2 * lx_matrix_csh[dim * (m1m + L) + (m2p + L)];
68 0 : lx += conj(b1) * a2 * lx_matrix_csh[dim * (m1p + L) + (m2m + L)];
69 0 : lx += conj(b1) * b2 * lx_matrix_csh[dim * (m1p + L) + (m2p + L)];
70 :
71 0 : double complex ly = 0.0 + 0.0 * I;
72 0 : ly += conj(a1) * a2 * ly_matrix_csh[dim * (m1m + L) + (m2m + L)];
73 0 : ly += conj(a1) * b2 * ly_matrix_csh[dim * (m1m + L) + (m2p + L)];
74 0 : ly += conj(b1) * a2 * ly_matrix_csh[dim * (m1p + L) + (m2m + L)];
75 0 : ly += conj(b1) * b2 * ly_matrix_csh[dim * (m1p + L) + (m2p + L)];
76 :
77 0 : double complex lz = 0.0 + 0.0 * I;
78 0 : lz += conj(a1) * a2 * lz_matrix_csh[dim * (m1m + L) + (m2m + L)];
79 0 : lz += conj(a1) * b2 * lz_matrix_csh[dim * (m1m + L) + (m2p + L)];
80 0 : lz += conj(b1) * a2 * lz_matrix_csh[dim * (m1p + L) + (m2m + L)];
81 0 : lz += conj(b1) * b2 * lz_matrix_csh[dim * (m1p + L) + (m2p + L)];
82 :
83 0 : lx_matrix[(m1 + L) * dim + (m2 + L)] = cimag(lx);
84 0 : ly_matrix[(m1 + L) * dim + (m2 + L)] = -creal(ly);
85 0 : lz_matrix[(m1 + L) * dim + (m2 + L)] = cimag(lz);
86 : }
87 : }
88 :
89 0 : free(lx_matrix_csh);
90 0 : free(ly_matrix_csh);
91 0 : free(lz_matrix_csh);
92 0 : }
93 :
94 : /**
95 : * Constructs matrices of the Lx, Ly, Lz operators in the basis of
96 : * complex spherical harmonics (csh) |Y_lm> for angular momentum l=L.
97 : * Matrices of size (2*L+1) x (2*L+1) must be pre-allocated.
98 : */
99 0 : void libgrpp_construct_angular_momentum_matrices_csh(int L, double *lx_matrix,
100 : double *ly_matrix,
101 : double *lz_matrix) {
102 0 : int dim = 2 * L + 1;
103 :
104 : // set all matrices to zero
105 0 : memset(lx_matrix, 0, dim * dim * sizeof(double));
106 0 : memset(ly_matrix, 0, dim * dim * sizeof(double));
107 0 : memset(lz_matrix, 0, dim * dim * sizeof(double));
108 :
109 0 : for (int m1 = -L; m1 <= L; m1++) {
110 0 : for (int m2 = -L; m2 <= L; m2++) {
111 :
112 0 : double lz = m2 * (m1 == m2);
113 0 : double lp = sqrt((L - m2) * (L + m2 + 1)) * (m1 == m2 + 1); // L+
114 0 : double lm = sqrt((L + m2) * (L - m2 + 1)) * (m1 == m2 - 1); // L-
115 0 : double lx = 0.5 * (lp + lm);
116 0 : double ly = 0.5 * (lp - lm);
117 :
118 0 : lx_matrix[(m1 + L) * dim + (m2 + L)] = lx;
119 0 : ly_matrix[(m1 + L) * dim + (m2 + L)] = ly;
120 0 : lz_matrix[(m1 + L) * dim + (m2 + L)] = lz;
121 : }
122 : }
123 0 : }
124 :
125 : /**
126 : * Real spherical harmonic S_{l,m} can be represented as a linear combination
127 : * of two complex spherical harmonics:
128 : * S_{l,m} = a * Y_{l,-m} + b * Y_{l,m}
129 : * (except for the case m=0, where S_{l,0} = Y_{l,0})
130 : *
131 : * coefficients can be found elsewhere, see, for example,
132 : * https://en.wikipedia.org/wiki/Table_of_spherical_harmonics
133 : */
134 0 : static void get_transformation_coeffs_csh_to_rsh(int m, double complex *a,
135 : double complex *b) {
136 0 : if (m == 0) {
137 0 : *a = 0.5;
138 0 : *b = 0.5;
139 0 : } else if (m < 0) {
140 0 : *a = +1.0 * I / sqrt(2);
141 0 : *b = -1.0 * I / sqrt(2) * pow(-1, abs(m));
142 : } else { // m > 0
143 0 : *a = +1.0 / sqrt(2);
144 0 : *b = +1.0 / sqrt(2) * pow(-1, m);
145 : }
146 0 : }
|
