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: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include <stdbool.h>
9 :
10 : #if defined(__CUDACC__) || defined(__HIPCC__)
11 : #define GRID_DEVICE __device__
12 : #else
13 : #define GRID_DEVICE
14 : #endif
15 :
16 : /*******************************************************************************
17 : * \brief Returns matrix element cab[idx(b)][idx(a)].
18 : * This function has to be implemented by the importing compilation unit.
19 : * A simple implementation is just: returns cab[idx(b) * n1 + idx(a)];
20 : * \author Ole Schuett
21 : ******************************************************************************/
22 : GRID_DEVICE static inline double cab_get(const cab_store *cab, const orbital a,
23 : const orbital b);
24 :
25 : /*******************************************************************************
26 : * \brief Returns i'th component of force on atom a for compute_tau=false.
27 : * \author Ole Schuett
28 : ******************************************************************************/
29 : GRID_DEVICE static inline double
30 467882250 : get_force_a_normal(const orbital a, const orbital b, const int i,
31 : const double zeta, const cab_store *cab) {
32 467882250 : const double aip1 = cab_get(cab, up(i, a), b);
33 467882250 : const double aim1 = cab_get(cab, down(i, a), b);
34 467882250 : return 2.0 * zeta * aip1 - a.l[i] * aim1;
35 : }
36 :
37 : /*******************************************************************************
38 : * \brief Returns i'th component of force on atom a.
39 : * \author Ole Schuett
40 : ******************************************************************************/
41 : GRID_DEVICE static inline double
42 431261952 : get_force_a(const orbital a, const orbital b, const int i, const double zeta,
43 : const double zetb, const cab_store *cab, const bool compute_tau) {
44 431261952 : if (!compute_tau) {
45 427932834 : return get_force_a_normal(a, b, i, zeta, cab);
46 : } else {
47 : double force = 0.0;
48 13316472 : for (int k = 0; k < 3; k++) {
49 9987354 : force += 0.5 * a.l[k] * b.l[k] *
50 9987354 : get_force_a_normal(down(k, a), down(k, b), i, zeta, cab);
51 9987354 : force -= zeta * b.l[k] *
52 9987354 : get_force_a_normal(up(k, a), down(k, b), i, zeta, cab);
53 9987354 : force -= a.l[k] * zetb *
54 9987354 : get_force_a_normal(down(k, a), up(k, b), i, zeta, cab);
55 9987354 : force += 2.0 * zeta * zetb *
56 9987354 : get_force_a_normal(up(k, a), up(k, b), i, zeta, cab);
57 : }
58 : return force;
59 : }
60 : }
61 :
62 : /*******************************************************************************
63 : * \brief Returns i'th component of force on atom b for compute_tau=false.
64 : * \author Ole Schuett
65 : ******************************************************************************/
66 : GRID_DEVICE static inline double
67 467853894 : get_force_b_normal(const orbital a, const orbital b, const int i,
68 : const double zetb, const double rab[3],
69 : const cab_store *cab) {
70 467853894 : const double axpm0 = cab_get(cab, a, b);
71 467853894 : const double aip1 = cab_get(cab, up(i, a), b);
72 467853894 : const double bim1 = cab_get(cab, a, down(i, b));
73 467853894 : return 2.0 * zetb * (aip1 - rab[i] * axpm0) - b.l[i] * bim1;
74 : }
75 :
76 : /*******************************************************************************
77 : * \brief Returns i'th component of force on atom b.
78 : * \author Ole Schuett
79 : ******************************************************************************/
80 : GRID_DEVICE static inline double
81 431233596 : get_force_b(const orbital a, const orbital b, const int i, const double zeta,
82 : const double zetb, const double rab[3], const cab_store *cab,
83 : const bool compute_tau) {
84 431233596 : if (!compute_tau) {
85 427904478 : return get_force_b_normal(a, b, i, zetb, rab, cab);
86 : } else {
87 : double force = 0.0;
88 13316472 : for (int k = 0; k < 3; k++) {
89 9987354 : force += 0.5 * a.l[k] * b.l[k] *
90 9987354 : get_force_b_normal(down(k, a), down(k, b), i, zetb, rab, cab);
91 9987354 : force -= zeta * b.l[k] *
92 9987354 : get_force_b_normal(up(k, a), down(k, b), i, zetb, rab, cab);
93 9987354 : force -= a.l[k] * zetb *
94 9987354 : get_force_b_normal(down(k, a), up(k, b), i, zetb, rab, cab);
95 9987354 : force += 2.0 * zeta * zetb *
96 9987354 : get_force_b_normal(up(k, a), up(k, b), i, zetb, rab, cab);
97 : }
98 : return force;
99 : }
100 : }
101 :
102 : /*******************************************************************************
103 : * \brief Returns element i,j of virial on atom a for compute_tau=false.
104 : * \author Ole Schuett
105 : ******************************************************************************/
106 : GRID_DEVICE static inline double
107 264337245 : get_virial_a_normal(const orbital a, const orbital b, const int i, const int j,
108 : const double zeta, const cab_store *cab) {
109 264337245 : return 2.0 * zeta * cab_get(cab, up(i, up(j, a)), b) -
110 264337245 : a.l[j] * cab_get(cab, up(i, down(j, a)), b);
111 : }
112 :
113 : /*******************************************************************************
114 : * \brief Returns element i,j of virial on atom a.
115 : * \author Ole Schuett
116 : ******************************************************************************/
117 : GRID_DEVICE static inline double
118 199520757 : get_virial_a(const orbital a, const orbital b, const int i, const int j,
119 : const double zeta, const double zetb, const cab_store *cab,
120 : const bool compute_tau) {
121 :
122 199520757 : if (!compute_tau) {
123 193628349 : return get_virial_a_normal(a, b, i, j, zeta, cab);
124 : } else {
125 : double virial = 0.0;
126 23569632 : for (int k = 0; k < 3; k++) {
127 17677224 : virial += 0.5 * a.l[k] * b.l[k] *
128 17677224 : get_virial_a_normal(down(k, a), down(k, b), i, j, zeta, cab);
129 17677224 : virial -= zeta * b.l[k] *
130 17677224 : get_virial_a_normal(up(k, a), down(k, b), i, j, zeta, cab);
131 17677224 : virial -= a.l[k] * zetb *
132 17677224 : get_virial_a_normal(down(k, a), up(k, b), i, j, zeta, cab);
133 17677224 : virial += 2.0 * zeta * zetb *
134 17677224 : get_virial_a_normal(up(k, a), up(k, b), i, j, zeta, cab);
135 : }
136 : return virial;
137 : }
138 : }
139 :
140 : /*******************************************************************************
141 : * \brief Returns element i,j of virial on atom b for compute_tau=false.
142 : * \author Ole Schuett
143 : ******************************************************************************/
144 : GRID_DEVICE static inline double
145 264334005 : get_virial_b_normal(const orbital a, const orbital b, const int i, const int j,
146 : const double zetb, const double rab[3],
147 : const cab_store *cab) {
148 :
149 264334005 : return 2.0 * zetb *
150 264334005 : (cab_get(cab, up(i, up(j, a)), b) -
151 264334005 : cab_get(cab, up(i, a), b) * rab[j] -
152 264334005 : cab_get(cab, up(j, a), b) * rab[i] +
153 264334005 : cab_get(cab, a, b) * rab[j] * rab[i]) -
154 264334005 : b.l[j] * cab_get(cab, a, up(i, down(j, b)));
155 : }
156 :
157 : /*******************************************************************************
158 : * \brief Returns element i,j of virial on atom b.
159 : * \author Ole Schuett
160 : ******************************************************************************/
161 : GRID_DEVICE static inline double
162 199517517 : get_virial_b(const orbital a, const orbital b, const int i, const int j,
163 : const double zeta, const double zetb, const double rab[3],
164 : const cab_store *cab, const bool compute_tau) {
165 :
166 199517517 : if (!compute_tau) {
167 193625109 : return get_virial_b_normal(a, b, i, j, zetb, rab, cab);
168 : } else {
169 : double virial = 0.0;
170 23569632 : for (int k = 0; k < 3; k++) {
171 17677224 : virial +=
172 17677224 : 0.5 * a.l[k] * b.l[k] *
173 17677224 : get_virial_b_normal(down(k, a), down(k, b), i, j, zetb, rab, cab);
174 17677224 : virial -= zeta * b.l[k] *
175 17677224 : get_virial_b_normal(up(k, a), down(k, b), i, j, zetb, rab, cab);
176 17677224 : virial -= a.l[k] * zetb *
177 17677224 : get_virial_b_normal(down(k, a), up(k, b), i, j, zetb, rab, cab);
178 17677224 : virial += 2.0 * zeta * zetb *
179 17677224 : get_virial_b_normal(up(k, a), up(k, b), i, j, zetb, rab, cab);
180 : }
181 : return virial;
182 : }
183 : }
184 :
185 : /*******************************************************************************
186 : * \brief Returns element i,j of hab matrix.
187 : * \author Ole Schuett
188 : ******************************************************************************/
189 1053102520 : GRID_DEVICE static inline double get_hab(const orbital a, const orbital b,
190 : const double zeta, const double zetb,
191 : const cab_store *cab,
192 : const bool compute_tau) {
193 1053102520 : if (!compute_tau) {
194 1042171833 : return cab_get(cab, a, b);
195 : } else {
196 : double hab = 0.0;
197 43722748 : for (int k = 0; k < 3; k++) {
198 32792061 : hab += 0.5 * a.l[k] * b.l[k] * cab_get(cab, down(k, a), down(k, b));
199 32792061 : hab -= zeta * b.l[k] * cab_get(cab, up(k, a), down(k, b));
200 32792061 : hab -= a.l[k] * zetb * cab_get(cab, down(k, a), up(k, b));
201 32792061 : hab += 2.0 * zeta * zetb * cab_get(cab, up(k, a), up(k, b));
202 : }
203 : return hab;
204 : }
205 : }
206 :
207 : /*******************************************************************************
208 : * \brief Differences in angular momentum.
209 : * \author Ole Schuett
210 : ******************************************************************************/
211 : typedef struct {
212 : int la_max_diff;
213 : int la_min_diff;
214 : int lb_max_diff;
215 : int lb_min_diff;
216 : } process_ldiffs;
217 :
218 : /*******************************************************************************
219 : * \brief Returns difference in angular momentum range for given flags.
220 : * \author Ole Schuett
221 : ******************************************************************************/
222 90811653 : static process_ldiffs process_get_ldiffs(bool calculate_forces,
223 : bool calculate_virial,
224 : bool compute_tau) {
225 90811653 : process_ldiffs ldiffs;
226 :
227 90811653 : ldiffs.la_max_diff = 0;
228 90811653 : ldiffs.lb_max_diff = 0;
229 90811653 : ldiffs.la_min_diff = 0;
230 90811653 : ldiffs.lb_min_diff = 0;
231 :
232 90811653 : if (calculate_forces || calculate_virial) {
233 18660552 : ldiffs.la_max_diff += 1; // for deriv. of gaussian, unimportant which one
234 18660552 : ldiffs.la_min_diff -= 1;
235 18660552 : ldiffs.lb_min_diff -= 1;
236 : }
237 :
238 18660552 : if (calculate_virial) {
239 3024467 : ldiffs.la_max_diff += 1;
240 3024467 : ldiffs.lb_max_diff += 1;
241 : }
242 :
243 90811653 : if (compute_tau) {
244 1353452 : ldiffs.la_max_diff += 1;
245 1353452 : ldiffs.lb_max_diff += 1;
246 1353452 : ldiffs.la_min_diff -= 1;
247 1353452 : ldiffs.lb_min_diff -= 1;
248 : }
249 :
250 90811653 : return ldiffs;
251 : }
252 :
253 : // EOF
|
