Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2023 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 : * \author Ole Schuett
19 : ******************************************************************************/
20 4977913052 : GRID_DEVICE static inline double get_term(const orbital a, const orbital b,
21 : const int n, const double *cab) {
22 4977913052 : return cab[idx(b) * n + idx(a)];
23 : }
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 449265999 : get_force_a_normal(const orbital a, const orbital b, const int i,
31 : const double zeta, const int n, const double *cab) {
32 449265999 : const double aip1 = get_term(up(i, a), b, n, cab);
33 449265999 : const double aim1 = get_term(down(i, a), b, n, cab);
34 449265999 : 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 413988141 : GRID_DEVICE static inline double get_force_a(const orbital a, const orbital b,
42 : const int i, const double zeta,
43 : const double zetb, const int n,
44 : const double *cab,
45 : const bool compute_tau) {
46 413988141 : if (!compute_tau) {
47 410781063 : return get_force_a_normal(a, b, i, zeta, n, cab);
48 : } else {
49 : double force = 0.0;
50 12828312 : for (int k = 0; k < 3; k++) {
51 9621234 : force += 0.5 * a.l[k] * b.l[k] *
52 9621234 : get_force_a_normal(down(k, a), down(k, b), i, zeta, n, cab);
53 9621234 : force -= zeta * b.l[k] *
54 9621234 : get_force_a_normal(up(k, a), down(k, b), i, zeta, n, cab);
55 9621234 : force -= a.l[k] * zetb *
56 9621234 : get_force_a_normal(down(k, a), up(k, b), i, zeta, n, cab);
57 9621234 : force += 2.0 * zeta * zetb *
58 9621234 : get_force_a_normal(up(k, a), up(k, b), i, zeta, n, cab);
59 : }
60 : return force;
61 : }
62 : }
63 :
64 : /*******************************************************************************
65 : * \brief Returns i'th component of force on atom b for compute_tau=false.
66 : * \author Ole Schuett
67 : ******************************************************************************/
68 : GRID_DEVICE static inline double
69 449250267 : get_force_b_normal(const orbital a, const orbital b, const int i,
70 : const double zetb, const double rab[3], const int n,
71 : const double *cab) {
72 449250267 : const double axpm0 = get_term(a, b, n, cab);
73 449250267 : const double aip1 = get_term(up(i, a), b, n, cab);
74 449250267 : const double bim1 = get_term(a, down(i, b), n, cab);
75 449250267 : return 2.0 * zetb * (aip1 - rab[i] * axpm0) - b.l[i] * bim1;
76 : }
77 :
78 : /*******************************************************************************
79 : * \brief Returns i'th component of force on atom b.
80 : * \author Ole Schuett
81 : ******************************************************************************/
82 : GRID_DEVICE static inline double
83 413972409 : get_force_b(const orbital a, const orbital b, const int i, const double zeta,
84 : const double zetb, const double rab[3], const int n,
85 : const double *cab, const bool compute_tau) {
86 413972409 : if (!compute_tau) {
87 410765331 : return get_force_b_normal(a, b, i, zetb, rab, n, cab);
88 : } else {
89 : double force = 0.0;
90 12828312 : for (int k = 0; k < 3; k++) {
91 9621234 : force += 0.5 * a.l[k] * b.l[k] *
92 9621234 : get_force_b_normal(down(k, a), down(k, b), i, zetb, rab, n, cab);
93 9621234 : force -= zeta * b.l[k] *
94 9621234 : get_force_b_normal(up(k, a), down(k, b), i, zetb, rab, n, cab);
95 9621234 : force -= a.l[k] * zetb *
96 9621234 : get_force_b_normal(down(k, a), up(k, b), i, zetb, rab, n, cab);
97 9621234 : force += 2.0 * zeta * zetb *
98 9621234 : get_force_b_normal(up(k, a), up(k, b), i, zetb, rab, n, cab);
99 : }
100 : return force;
101 : }
102 : }
103 :
104 : /*******************************************************************************
105 : * \brief Returns element i,j of virial on atom a for compute_tau=false.
106 : * \author Ole Schuett
107 : ******************************************************************************/
108 : GRID_DEVICE static inline double
109 249402600 : get_virial_a_normal(const orbital a, const orbital b, const int i, const int j,
110 : const double zeta, const int n, const double *cab) {
111 249402600 : return 2.0 * zeta * get_term(up(i, up(j, a)), b, n, cab) -
112 249402600 : a.l[j] * get_term(up(i, down(j, a)), b, n, cab);
113 : }
114 :
115 : /*******************************************************************************
116 : * \brief Returns element i,j of virial on atom a.
117 : * \author Ole Schuett
118 : ******************************************************************************/
119 : GRID_DEVICE static inline double
120 186596208 : get_virial_a(const orbital a, const orbital b, const int i, const int j,
121 : const double zeta, const double zetb, const int n,
122 : const double *cab, const bool compute_tau) {
123 :
124 186596208 : if (!compute_tau) {
125 180886536 : return get_virial_a_normal(a, b, i, j, zeta, n, cab);
126 : } else {
127 : double virial = 0.0;
128 22838688 : for (int k = 0; k < 3; k++) {
129 17129016 : virial += 0.5 * a.l[k] * b.l[k] *
130 17129016 : get_virial_a_normal(down(k, a), down(k, b), i, j, zeta, n, cab);
131 17129016 : virial -= zeta * b.l[k] *
132 17129016 : get_virial_a_normal(up(k, a), down(k, b), i, j, zeta, n, cab);
133 17129016 : virial -= a.l[k] * zetb *
134 17129016 : get_virial_a_normal(down(k, a), up(k, b), i, j, zeta, n, cab);
135 17129016 : virial += 2.0 * zeta * zetb *
136 17129016 : get_virial_a_normal(up(k, a), up(k, b), i, j, zeta, n, cab);
137 : }
138 : return virial;
139 : }
140 : }
141 :
142 : /*******************************************************************************
143 : * \brief Returns element i,j of virial on atom b for compute_tau=false.
144 : * \author Ole Schuett
145 : ******************************************************************************/
146 : GRID_DEVICE static inline double
147 249399720 : get_virial_b_normal(const orbital a, const orbital b, const int i, const int j,
148 : const double zetb, const double rab[3], const int n,
149 : const double *cab) {
150 :
151 249399720 : return 2.0 * zetb *
152 249399720 : (get_term(up(i, up(j, a)), b, n, cab) -
153 249399720 : get_term(up(i, a), b, n, cab) * rab[j] -
154 249399720 : get_term(up(j, a), b, n, cab) * rab[i] +
155 249399720 : get_term(a, b, n, cab) * rab[j] * rab[i]) -
156 249399720 : b.l[j] * get_term(a, up(i, down(j, b)), n, cab);
157 : }
158 :
159 : /*******************************************************************************
160 : * \brief Returns element i,j of virial on atom b.
161 : * \author Ole Schuett
162 : ******************************************************************************/
163 : GRID_DEVICE static inline double
164 186593328 : get_virial_b(const orbital a, const orbital b, const int i, const int j,
165 : const double zeta, const double zetb, const double rab[3],
166 : const int n, const double *cab, const bool compute_tau) {
167 :
168 186593328 : if (!compute_tau) {
169 180883656 : return get_virial_b_normal(a, b, i, j, zetb, rab, n, cab);
170 : } else {
171 : double virial = 0.0;
172 22838688 : for (int k = 0; k < 3; k++) {
173 17129016 : virial +=
174 17129016 : 0.5 * a.l[k] * b.l[k] *
175 17129016 : get_virial_b_normal(down(k, a), down(k, b), i, j, zetb, rab, n, cab);
176 17129016 : virial -=
177 17129016 : zeta * b.l[k] *
178 17129016 : get_virial_b_normal(up(k, a), down(k, b), i, j, zetb, rab, n, cab);
179 17129016 : virial -=
180 17129016 : a.l[k] * zetb *
181 17129016 : get_virial_b_normal(down(k, a), up(k, b), i, j, zetb, rab, n, cab);
182 17129016 : virial +=
183 17129016 : 2.0 * zeta * zetb *
184 17129016 : get_virial_b_normal(up(k, a), up(k, b), i, j, zetb, rab, n, cab);
185 : }
186 : return virial;
187 : }
188 : }
189 :
190 : /*******************************************************************************
191 : * \brief Returns element i,j of hab matrix.
192 : * \author Ole Schuett
193 : ******************************************************************************/
194 873879761 : GRID_DEVICE static inline double get_hab(const orbital a, const orbital b,
195 : const double zeta, const double zetb,
196 : const int n, const double *cab,
197 : const bool compute_tau) {
198 873879761 : if (!compute_tau) {
199 863702789 : return get_term(a, b, n, cab);
200 : } else {
201 : double hab = 0.0;
202 40707888 : for (int k = 0; k < 3; k++) {
203 30530916 : hab += 0.5 * a.l[k] * b.l[k] * get_term(down(k, a), down(k, b), n, cab);
204 30530916 : hab -= zeta * b.l[k] * get_term(up(k, a), down(k, b), n, cab);
205 30530916 : hab -= a.l[k] * zetb * get_term(down(k, a), up(k, b), n, cab);
206 30530916 : hab += 2.0 * zeta * zetb * get_term(up(k, a), up(k, b), n, cab);
207 : }
208 : return hab;
209 : }
210 : }
211 :
212 : /*******************************************************************************
213 : * \brief Differences in angular momentum.
214 : * \author Ole Schuett
215 : ******************************************************************************/
216 : typedef struct {
217 : int la_max_diff;
218 : int la_min_diff;
219 : int lb_max_diff;
220 : int lb_min_diff;
221 : } process_ldiffs;
222 :
223 : /*******************************************************************************
224 : * \brief Returns difference in angular momentum range for given flags.
225 : * \author Ole Schuett
226 : ******************************************************************************/
227 85258416 : static process_ldiffs process_get_ldiffs(bool calculate_forces,
228 : bool calculate_virial,
229 : bool compute_tau) {
230 85258416 : process_ldiffs ldiffs;
231 :
232 85258416 : ldiffs.la_max_diff = 0;
233 85258416 : ldiffs.lb_max_diff = 0;
234 85258416 : ldiffs.la_min_diff = 0;
235 85258416 : ldiffs.lb_min_diff = 0;
236 :
237 85258416 : if (calculate_forces || calculate_virial) {
238 17692668 : ldiffs.la_max_diff += 1; // for deriv. of gaussian, unimportant which one
239 17692668 : ldiffs.la_min_diff -= 1;
240 17692668 : ldiffs.lb_min_diff -= 1;
241 : }
242 :
243 85258416 : if (calculate_virial) {
244 2730085 : ldiffs.la_max_diff += 1;
245 2730085 : ldiffs.lb_max_diff += 1;
246 : }
247 :
248 85258416 : if (compute_tau) {
249 1230687 : ldiffs.la_max_diff += 1;
250 1230687 : ldiffs.lb_max_diff += 1;
251 1230687 : ldiffs.la_min_diff -= 1;
252 1230687 : ldiffs.lb_min_diff -= 1;
253 : }
254 :
255 85258416 : return ldiffs;
256 : }
257 :
258 : // EOF
|
