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 : #include <stdlib.h>
16 : #include <string.h>
17 :
18 : #include "grpp_overlap_gradient.h"
19 : #include "libgrpp.h"
20 :
21 : #include "grpp_diff_gaussian.h"
22 : #include "grpp_utils.h"
23 :
24 : static void overlap_gradient_diff_bra_contribution(libgrpp_shell_t *shell_A,
25 : libgrpp_shell_t *shell_B,
26 : double **grad,
27 : double factor);
28 :
29 : static void overlap_gradient_diff_bra_overlap_integrals(
30 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double **overlap_down,
31 : double **overlap_up, int *cart_size_down, int *cart_size_up);
32 :
33 : int nlm_to_linear(int *nlm);
34 :
35 : /**
36 : * Analytic calculation of gradients of overlap integrals for a given shell pair
37 : * with respect to the point 'point_3d'.
38 : */
39 0 : void libgrpp_overlap_integrals_gradient(libgrpp_shell_t *shell_A,
40 : libgrpp_shell_t *shell_B,
41 : double *point_3d, double **grad) {
42 0 : int cart_size_A = libgrpp_get_shell_size(shell_A);
43 0 : int cart_size_B = libgrpp_get_shell_size(shell_B);
44 0 : int buf_size = cart_size_A * cart_size_B;
45 :
46 : /*
47 : * initializations: set gradients to zero
48 : */
49 0 : for (int icoord = 0; icoord < 3; icoord++) {
50 0 : memset(grad[icoord], 0, sizeof(double) * buf_size);
51 : }
52 :
53 : /*
54 : * integrals are zero:
55 : * (1) for 1-center integrals <A|A> (due to the translational invariance)
56 : * (2) d<A|B> / dC = 0 (integral is constant for the given 'point_3d')
57 : */
58 0 : if (points_are_equal(shell_A->origin, shell_B->origin)) {
59 : return;
60 : }
61 :
62 : /*
63 : * construct gradients:
64 : * d<A|B>/dA = + < df/dA | B >
65 : * d<A|B>/dB = - < df/dA | B >
66 : *
67 : * note that due to the property of translational invariance,
68 : * d<A|B>/dB = - d<A|B>/dA
69 : */
70 0 : if (points_are_equal(shell_A->origin, point_3d)) {
71 0 : overlap_gradient_diff_bra_contribution(shell_A, shell_B, grad, +1.0);
72 : }
73 0 : if (points_are_equal(shell_B->origin, point_3d)) {
74 0 : overlap_gradient_diff_bra_contribution(shell_A, shell_B, grad, -1.0);
75 : }
76 : }
77 :
78 : /**
79 : * Calculates contribution to gradients arising from the < df/dA | g > term:
80 : *
81 : * grad += factor * < df/dA | g >
82 : *
83 : * (bra basis function is differentiated).
84 : */
85 0 : static void overlap_gradient_diff_bra_contribution(libgrpp_shell_t *shell_A,
86 : libgrpp_shell_t *shell_B,
87 : double **grad,
88 : double factor) {
89 : /*
90 : * calculate overlap integrals < df/dA | B >
91 : */
92 0 : double *overlap_down = NULL;
93 0 : double *overlap_up = NULL;
94 0 : int cart_size_down = 0;
95 0 : int cart_size_up = 0;
96 :
97 0 : overlap_gradient_diff_bra_overlap_integrals(shell_A, shell_B, &overlap_down,
98 : &overlap_up, &cart_size_down,
99 : &cart_size_up);
100 :
101 : /*
102 : * construct contributions to gradients:
103 : * d<A|B>/dA += < df/dA | B >
104 : */
105 0 : for (int icoord = 0; icoord < 3; icoord++) {
106 0 : for (int i = 0; i < shell_A->cart_size; i++) {
107 0 : for (int j = 0; j < shell_B->cart_size; j++) {
108 :
109 0 : int *bra_nlm = shell_A->cart_list + 3 * i;
110 0 : int *ket_nlm = shell_B->cart_list + 3 * j;
111 0 : int index = i * shell_B->cart_size + j;
112 :
113 : /*
114 : * contribution from the L-1 gaussian
115 : */
116 0 : if (shell_A->L > 0) {
117 0 : bra_nlm[icoord] -= 1;
118 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
119 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
120 0 : bra_nlm[icoord] += 1;
121 :
122 0 : grad[icoord][index] -=
123 0 : factor * bra_nlm[icoord] *
124 0 : overlap_down[shell_B->cart_size * bra_index + ket_index];
125 : }
126 :
127 : /*
128 : * contribution from the L+1 gaussian
129 : */
130 0 : bra_nlm[icoord] += 1;
131 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
132 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
133 0 : bra_nlm[icoord] -= 1;
134 :
135 0 : grad[icoord][index] +=
136 0 : factor * overlap_up[shell_B->cart_size * bra_index + ket_index];
137 : }
138 : }
139 : }
140 :
141 0 : if (overlap_down) {
142 0 : free(overlap_down);
143 : }
144 0 : free(overlap_up);
145 0 : }
146 :
147 : /**
148 : * To assemble the contribution < df/dA | g > to gradients, one have to
149 : * differentiate Gaussian function. Such a differentiation yields two Gaussians
150 : * with angular momenta L-1 ("down") and L+1 ("up"): dG/dA -> G(L-1) and G(L+1)
151 : *
152 : * This function constructs overlap matrices with these "downgraded" and
153 : * "upgraded" Gaussian functions: < G(L-1) | G' > and < G(L+1) | G' >
154 : *
155 : */
156 0 : static void overlap_gradient_diff_bra_overlap_integrals(
157 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double **overlap_down,
158 : double **overlap_up, int *cart_size_down, int *cart_size_up) {
159 : /*
160 : * differentiation of contracted Gaussian functions
161 : */
162 0 : libgrpp_shell_t *shell_A_down = NULL;
163 0 : libgrpp_shell_t *shell_A_up = NULL;
164 0 : libgrpp_differentiate_shell(shell_A, &shell_A_down, &shell_A_up);
165 :
166 0 : *cart_size_down = 0;
167 0 : if (shell_A_down != NULL) {
168 0 : *cart_size_down = shell_A_down->cart_size;
169 : }
170 0 : *cart_size_up = shell_A_up->cart_size;
171 :
172 : /*
173 : * overlap matrix:
174 : * < L-1 | L>
175 : */
176 0 : if (shell_A_down != NULL) {
177 0 : *overlap_down = (double *)calloc(
178 0 : shell_A_down->cart_size * shell_B->cart_size, sizeof(double));
179 0 : libgrpp_overlap_integrals(shell_A_down, shell_B, *overlap_down);
180 : } else {
181 0 : *overlap_down = NULL;
182 : }
183 :
184 : /*
185 : * overlap matrix:
186 : * < L+1 | L>
187 : */
188 0 : *overlap_up = (double *)calloc(shell_A_up->cart_size * shell_B->cart_size,
189 : sizeof(double));
190 0 : libgrpp_overlap_integrals(shell_A_up, shell_B, *overlap_up);
191 :
192 : /*
193 : * clean up
194 : */
195 0 : if (shell_A_down) {
196 0 : libgrpp_delete_shell(shell_A_down);
197 : }
198 0 : libgrpp_delete_shell(shell_A_up);
199 0 : }
200 :
201 : /**
202 : * calculates sequential ("linear") index of the (n,l,m) primitive in the
203 : * cartesian shell
204 : */
205 0 : int libgrpp_nlm_to_linear(int *nlm) {
206 0 : int n = nlm[0];
207 0 : int l = nlm[1];
208 0 : int m = nlm[2];
209 :
210 0 : int L = n + l + m;
211 0 : int cart_size = (L + 1) * (L + 2) / 2;
212 0 : int *cart_list = libgrpp_generate_shell_cartesians(L);
213 :
214 0 : int index = 0;
215 0 : for (index = 0; index < cart_size; index++) {
216 0 : if (cart_list[3 * index + 0] == n && cart_list[3 * index + 1] == l &&
217 0 : cart_list[3 * index + 2] == m) {
218 : break;
219 : }
220 : }
221 :
222 0 : free(cart_list);
223 :
224 0 : return index;
225 : }
|
