Psi4
dft_integrators.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2018 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This file is part of Psi4.
12  *
13  * Psi4 is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser General Public License as published by
15  * the Free Software Foundation, version 3.
16  *
17  * Psi4 is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License along
23  * with Psi4; if not, write to the Free Software Foundation, Inc.,
24  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25  *
26  * @END LICENSE
27  */
28 
29 #ifndef LIBFOCK_DFT_INTEGRATORS
30 #define LIBFOCK_DFT_INTEGRATORS
31 
32 #include "cubature.h"
33 #include "points.h"
34 
37 #include "psi4/libqt/qt.h"
38 #include "psi4/psi4-dec.h"
39 
40 #include "psi4/libmints/basisset.h"
41 #include "psi4/libmints/matrix.h"
42 #include "psi4/libmints/vector.h"
43 
44 #include <cstdlib>
45 
46 #ifdef _OPENMP
47 #include <omp.h>
48 #endif
49 
50 namespace psi {
51 
52 namespace dft_integrators {
53 
54 inline std::vector<double> rks_quadrature_integrate(std::shared_ptr<BlockOPoints> block,
55  std::shared_ptr<SuperFunctional> fworker,
56  std::shared_ptr<PointFunctions> pworker) {
57  // Block data
58  int npoints = block->npoints();
59  double* x = block->x();
60  double* y = block->y();
61  double* z = block->z();
62  double* w = block->w();
63 
64  // Superfunctional data
65  double* zk = fworker->value("V")->pointer();
66  double* QTp = fworker->value("Q_TMP")->pointer();
67 
68  // Points data
69  double* rho_a = pworker->point_value("RHO_A")->pointer();
70 
71  // Build quadrature
72  std::vector<double> ret(5);
73  ret[0] = C_DDOT(npoints, w, 1, zk, 1);
74  for (int P = 0; P < npoints; P++) {
75  QTp[P] = w[P] * rho_a[P];
76  }
77  ret[1] = C_DDOT(npoints, w, 1, rho_a, 1);
78  ret[2] = C_DDOT(npoints, QTp, 1, x, 1);
79  ret[3] = C_DDOT(npoints, QTp, 1, y, 1);
80  ret[4] = C_DDOT(npoints, QTp, 1, z, 1);
81 
82  return ret;
83 }
84 
85 inline void rks_integrator(std::shared_ptr<BlockOPoints> block, std::shared_ptr<SuperFunctional> fworker,
86  std::shared_ptr<PointFunctions> pworker, SharedMatrix V, int ansatz = -1) {
87  ansatz = (ansatz == -1 ? fworker->ansatz() : ansatz);
88  // printf("Ansatz %d\n", ansatz);
89 
90  // Block data
91  const std::vector<int>& function_map = block->functions_local_to_global();
92  int nlocal = function_map.size();
93  int npoints = block->npoints();
94  double* w = block->w();
95 
96  // Scratch is updated
97  double** Tp = pworker->scratch()[0]->pointer();
98 
99  // Points data
100  double** phi = pworker->basis_value("PHI")->pointer();
101  double* rho_a = pworker->point_value("RHO_A")->pointer();
102  size_t coll_funcs = pworker->basis_value("PHI")->ncol();
103 
104  // V2 Temporary
105  int max_functions = V->ncol();
106  double** V2p = V->pointer();
107 
108  // => LSDA contribution (symmetrized) <= //
109  double* v_rho_a = fworker->value("V_RHO_A")->pointer();
110  for (int P = 0; P < npoints; P++) {
111  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
112  C_DAXPY(nlocal, 0.5 * v_rho_a[P] * w[P], phi[P], 1, Tp[P], 1);
113  }
114  // parallel_timer_off("LSDA Phi_tmp", rank);
115 
116  // => GGA contribution (symmetrized) <= //
117  if (ansatz >= 1) {
118  // parallel_timer_on("GGA Phi_tmp", rank);
119  double** phix = pworker->basis_value("PHI_X")->pointer();
120  double** phiy = pworker->basis_value("PHI_Y")->pointer();
121  double** phiz = pworker->basis_value("PHI_Z")->pointer();
122  double* rho_ax = pworker->point_value("RHO_AX")->pointer();
123  double* rho_ay = pworker->point_value("RHO_AY")->pointer();
124  double* rho_az = pworker->point_value("RHO_AZ")->pointer();
125  double* v_sigma_aa = fworker->value("V_GAMMA_AA")->pointer();
126 
127  for (int P = 0; P < npoints; P++) {
128  C_DAXPY(nlocal, w[P] * (2.0 * v_sigma_aa[P] * rho_ax[P]), phix[P], 1, Tp[P], 1);
129  C_DAXPY(nlocal, w[P] * (2.0 * v_sigma_aa[P] * rho_ay[P]), phiy[P], 1, Tp[P], 1);
130  C_DAXPY(nlocal, w[P] * (2.0 * v_sigma_aa[P] * rho_az[P]), phiz[P], 1, Tp[P], 1);
131  }
132  // parallel_timer_off("GGA Phi_tmp", rank);
133  }
134 
135  // Collect V terms
136  C_DGEMM('T', 'N', nlocal, nlocal, npoints, 1.0, phi[0], coll_funcs, Tp[0], max_functions, 0.0, V2p[0],
137  max_functions);
138 
139  for (int m = 0; m < nlocal; m++) {
140  for (int n = 0; n <= m; n++) {
141  V2p[m][n] = V2p[n][m] = V2p[m][n] + V2p[n][m];
142  }
143  }
144 
145  // => Meta contribution <= //
146  if (ansatz >= 2) {
147  // parallel_timer_on("Meta", rank);
148  double** phix = pworker->basis_value("PHI_X")->pointer();
149  double** phiy = pworker->basis_value("PHI_Y")->pointer();
150  double** phiz = pworker->basis_value("PHI_Z")->pointer();
151  double* v_tau_a = fworker->value("V_TAU_A")->pointer();
152 
153  double** phi_w[3];
154  phi_w[0] = phix;
155  phi_w[1] = phiy;
156  phi_w[2] = phiz;
157 
158  for (int i = 0; i < 3; i++) {
159  double** phiw = phi_w[i];
160  for (int P = 0; P < npoints; P++) {
161  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
162  C_DAXPY(nlocal, v_tau_a[P] * w[P], phiw[P], 1, Tp[P], 1);
163  }
164  C_DGEMM('T', 'N', nlocal, nlocal, npoints, 1.0, phiw[0], coll_funcs, Tp[0], max_functions, 1.0, V2p[0],
165  max_functions);
166  }
167  // parallel_timer_off("Meta", rank);
168  }
169 }
170 
171 inline void rks_gradient_integrator(std::shared_ptr<BasisSet> primary, std::shared_ptr<BlockOPoints> block,
172  std::shared_ptr<SuperFunctional> fworker, std::shared_ptr<PointFunctions> pworker,
173  SharedMatrix G, SharedMatrix U, int ansatz = -1) {
174  ansatz = (ansatz == -1 ? fworker->ansatz() : ansatz);
175 
176  // Get scratch pointers
177  double** Gp = G->pointer();
178 
179  double** Up = U->pointer();
180  double** Tp = pworker->scratch()[0]->pointer();
181  double** Dp = pworker->D_scratch()[0]->pointer();
182 
183  // Fine for now, but not true once we start caching
184  int max_functions = U->ncol();
185 
186  // Get block data
187  int npoints = block->npoints();
188  double* w = block->w();
189  const std::vector<int>& function_map = block->functions_local_to_global();
190  int nlocal = function_map.size();
191 
192  // Get points data
193  double** phi = pworker->basis_value("PHI")->pointer();
194  double** phi_x = pworker->basis_value("PHI_X")->pointer();
195  double** phi_y = pworker->basis_value("PHI_Y")->pointer();
196  double** phi_z = pworker->basis_value("PHI_Z")->pointer();
197  double* rho_a = pworker->point_value("RHO_A")->pointer();
198  size_t coll_funcs = pworker->basis_value("PHI")->ncol();
199 
200  // => LSDA Contribution <= //
201  double* v_rho_a = fworker->value("V_RHO_A")->pointer();
202  for (int P = 0; P < npoints; P++) {
203  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
204  C_DAXPY(nlocal, -2.0 * w[P] * v_rho_a[P], phi[P], 1, Tp[P], 1);
205  }
206 
207  // => GGA Contribution (Term 1) <= //
208  if (fworker->is_gga()) {
209  double* rho_ax = pworker->point_value("RHO_AX")->pointer();
210  double* rho_ay = pworker->point_value("RHO_AY")->pointer();
211  double* rho_az = pworker->point_value("RHO_AZ")->pointer();
212  double* v_gamma_aa = fworker->value("V_GAMMA_AA")->pointer();
213 
214  for (int P = 0; P < npoints; P++) {
215  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_ax[P]), phi_x[P], 1, Tp[P], 1);
216  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_ay[P]), phi_y[P], 1, Tp[P], 1);
217  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_az[P]), phi_z[P], 1, Tp[P], 1);
218  }
219  }
220 
221  // => Synthesis <= //
222  C_DGEMM('N', 'N', npoints, nlocal, nlocal, 1.0, Tp[0], max_functions, Dp[0], max_functions, 0.0, Up[0],
223  max_functions);
224 
225  for (int ml = 0; ml < nlocal; ml++) {
226  int A = primary->function_to_center(function_map[ml]);
227  Gp[A][0] += C_DDOT(npoints, &Up[0][ml], max_functions, &phi_x[0][ml], coll_funcs);
228  Gp[A][1] += C_DDOT(npoints, &Up[0][ml], max_functions, &phi_y[0][ml], coll_funcs);
229  Gp[A][2] += C_DDOT(npoints, &Up[0][ml], max_functions, &phi_z[0][ml], coll_funcs);
230  }
231 
232  // => GGA Contribution (Term 2) <= //
233  if (fworker->is_gga()) {
234  double** phi_xx = pworker->basis_value("PHI_XX")->pointer();
235  double** phi_xy = pworker->basis_value("PHI_XY")->pointer();
236  double** phi_xz = pworker->basis_value("PHI_XZ")->pointer();
237  double** phi_yy = pworker->basis_value("PHI_YY")->pointer();
238  double** phi_yz = pworker->basis_value("PHI_YZ")->pointer();
239  double** phi_zz = pworker->basis_value("PHI_ZZ")->pointer();
240  double* rho_ax = pworker->point_value("RHO_AX")->pointer();
241  double* rho_ay = pworker->point_value("RHO_AY")->pointer();
242  double* rho_az = pworker->point_value("RHO_AZ")->pointer();
243  double* v_gamma_aa = fworker->value("V_GAMMA_AA")->pointer();
244 
245  C_DGEMM('N', 'N', npoints, nlocal, nlocal, 1.0, phi[0], coll_funcs, Dp[0], max_functions, 0.0, Up[0],
246  max_functions);
247 
248  // x
249  for (int P = 0; P < npoints; P++) {
250  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
251  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_ax[P]), Up[P], 1, Tp[P], 1);
252  }
253  for (int ml = 0; ml < nlocal; ml++) {
254  int A = primary->function_to_center(function_map[ml]);
255  Gp[A][0] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_xx[0][ml], coll_funcs);
256  Gp[A][1] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_xy[0][ml], coll_funcs);
257  Gp[A][2] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_xz[0][ml], coll_funcs);
258  }
259 
260  // y
261  for (int P = 0; P < npoints; P++) {
262  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
263  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_ay[P]), Up[P], 1, Tp[P], 1);
264  }
265  for (int ml = 0; ml < nlocal; ml++) {
266  int A = primary->function_to_center(function_map[ml]);
267  Gp[A][0] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_xy[0][ml], coll_funcs);
268  Gp[A][1] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_yy[0][ml], coll_funcs);
269  Gp[A][2] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_yz[0][ml], coll_funcs);
270  }
271 
272  // z
273  for (int P = 0; P < npoints; P++) {
274  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
275  C_DAXPY(nlocal, -2.0 * w[P] * (2.0 * v_gamma_aa[P] * rho_az[P]), Up[P], 1, Tp[P], 1);
276  }
277  for (int ml = 0; ml < nlocal; ml++) {
278  int A = primary->function_to_center(function_map[ml]);
279  Gp[A][0] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_xz[0][ml], coll_funcs);
280  Gp[A][1] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_yz[0][ml], coll_funcs);
281  Gp[A][2] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_zz[0][ml], coll_funcs);
282  }
283  }
284 
285  // => Meta Contribution <= //
286  if (fworker->is_meta()) {
287  double** phi_xx = pworker->basis_value("PHI_XX")->pointer();
288  double** phi_xy = pworker->basis_value("PHI_XY")->pointer();
289  double** phi_xz = pworker->basis_value("PHI_XZ")->pointer();
290  double** phi_yy = pworker->basis_value("PHI_YY")->pointer();
291  double** phi_yz = pworker->basis_value("PHI_YZ")->pointer();
292  double** phi_zz = pworker->basis_value("PHI_ZZ")->pointer();
293  double* v_tau_a = fworker->value("V_TAU_A")->pointer();
294 
295  double** phi_i[3];
296  phi_i[0] = phi_x;
297  phi_i[1] = phi_y;
298  phi_i[2] = phi_z;
299 
300  double** phi_ij[3][3];
301  phi_ij[0][0] = phi_xx;
302  phi_ij[0][1] = phi_xy;
303  phi_ij[0][2] = phi_xz;
304  phi_ij[1][0] = phi_xy;
305  phi_ij[1][1] = phi_yy;
306  phi_ij[1][2] = phi_yz;
307  phi_ij[2][0] = phi_xz;
308  phi_ij[2][1] = phi_yz;
309  phi_ij[2][2] = phi_zz;
310 
311  for (int i = 0; i < 3; i++) {
312  double*** phi_j = phi_ij[i];
313  C_DGEMM('N', 'N', npoints, nlocal, nlocal, 1.0, phi_i[i][0], coll_funcs, Dp[0], max_functions, 0.0, Up[0],
314  max_functions);
315  for (int P = 0; P < npoints; P++) {
316  std::fill(Tp[P], Tp[P] + nlocal, 0.0);
317  C_DAXPY(nlocal, -2.0 * w[P] * (v_tau_a[P]), Up[P], 1, Tp[P], 1);
318  }
319  for (int ml = 0; ml < nlocal; ml++) {
320  int A = primary->function_to_center(function_map[ml]);
321  Gp[A][0] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_j[0][0][ml], coll_funcs);
322  Gp[A][1] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_j[1][0][ml], coll_funcs);
323  Gp[A][2] += C_DDOT(npoints, &Tp[0][ml], max_functions, &phi_j[2][0][ml], coll_funcs);
324  }
325  }
326  }
327 }
328 
329 } // End dft_integrator namespace
330 } // End psi namespace
331 #endif
int * U
Definition: stringlist.cc:66
void PSI_API C_DAXPY(size_t length, double a, double *x, int inc_x, double *y, int inc_y)
Definition: blas_intfc.cc:121
Definition: pointgrp.h:104
void rks_gradient_integrator(std::shared_ptr< BasisSet > primary, std::shared_ptr< BlockOPoints > block, std::shared_ptr< SuperFunctional > fworker, std::shared_ptr< PointFunctions > pworker, SharedMatrix G, SharedMatrix U, int ansatz=-1)
Definition: dft_integrators.h:171
Header file for the Quantum Trio LibraryDavid Sherrill 1994.
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
double PSI_API C_DDOT(size_t length, double *x, int inc_x, double *y, int inc_y)
Definition: blas_intfc.cc:214
PSI_API void C_DGEMM(char transa, char transb, int m, int n, int k, double alpha, double *a, int lda, double *b, int ldb, double beta, double *c, int ldc)
Definition: blas_intfc23.cc:324
std::vector< double > rks_quadrature_integrate(std::shared_ptr< BlockOPoints > block, std::shared_ptr< SuperFunctional > fworker, std::shared_ptr< PointFunctions > pworker)
Definition: dft_integrators.h:54
void rks_integrator(std::shared_ptr< BlockOPoints > block, std::shared_ptr< SuperFunctional > fworker, std::shared_ptr< PointFunctions > pworker, SharedMatrix V, int ansatz=-1)
Definition: dft_integrators.h:85