1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2011
9 * Copyright (c) 2010-2011, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #include "emdcube.h"
18 #include "gto_fourier.h"
19 #include "../basis.h"
20
21 #ifdef _OPENMP
22 #include <omp.h>
23 #endif
24
emd_cube(const BasisSet & bas,const arma::cx_mat & P,const std::vector<double> & px_arr,const std::vector<double> & py_arr,const std::vector<double> & pz_arr)25 void emd_cube(const BasisSet & bas, const arma::cx_mat & P, const std::vector<double> & px_arr, const std::vector<double> & py_arr, const std::vector<double> & pz_arr) {
26 // Get expansions of functions
27 std::vector< std::vector<size_t> > idents;
28 std::vector< std::vector<GTO_Fourier> > fourier=fourier_expand(bas,idents);
29
30 // Open output file.
31 FILE *out=fopen("emdcube.dat","w");
32
33 // Compute the norm (assumes evenly spaced grid!)
34 double norm=0.0;
35
36 // Compute the momentum densities in batches, allowing
37 // parallellization.
38 #ifdef _OPENMP
39 // The number of points per batch
40 const size_t Nbatch_p=100*omp_get_max_threads();
41 #else
42 const size_t Nbatch_p=100;
43 #endif
44
45 // The total number of point is
46 const size_t N=px_arr.size()*py_arr.size()*pz_arr.size();
47 // The necessary amount of batches is
48 size_t Nbatch=N/Nbatch_p;
49 if(N%Nbatch_p!=0)
50 Nbatch++;
51
52 // The values of momentum in the batch
53 coords_t p[Nbatch_p];
54 // The values of EMD in the batch
55 double emd[Nbatch_p];
56
57 // Number of points to compute in the batch
58 size_t np;
59 // Total number of points computed
60 size_t ntot=0;
61 // Indices of x, y and z
62 size_t xind=0, yind=0, zind=0;
63
64 // Loop over batches.
65 for(size_t ib=0;ib<Nbatch;ib++) {
66 // Zero amount of points in current batch.
67 np=0;
68
69 // Form list of points to compute.
70 while(np<Nbatch_p && ntot+np<N) {
71 p[np].x=px_arr[xind];
72 p[np].y=py_arr[yind];
73 p[np].z=pz_arr[zind];
74
75 // Increment number of points
76 np++;
77
78 // Determine next point.
79 if(zind+1<pz_arr.size())
80 zind++;
81 else {
82 zind=0;
83
84 if(yind+1<py_arr.size())
85 yind++;
86 else {
87 yind=0;
88 xind++;
89 }
90 }
91 }
92
93 // Begin parallel section
94 #ifdef _OPENMP
95 #pragma omp parallel
96 #endif
97 {
98
99 #ifdef _OPENMP
100 #pragma omp for
101 #endif
102 // Loop over the points in the batch
103 for(size_t ip=0;ip<np;ip++)
104 // Evaluate EMD
105 emd[ip]=eval_emd(bas,P,fourier,idents,p[ip].x,p[ip].y,p[ip].z);
106 } // end parallel region
107
108 // Save computed value of EMD and increment norm
109 for(size_t ip=0;ip<np;ip++) {
110 fprintf(out,"%e\t%e\t%e\t%e\n",p[ip].x,p[ip].y,p[ip].z,emd[ip]);
111 norm+=emd[ip];
112 }
113
114 // Increment number of computed points
115 ntot+=np;
116 }
117 // Close output file.
118 fclose(out);
119
120 // Plug in the spacing in the integral
121 double dx=(px_arr[px_arr.size()-1]-px_arr[0])/px_arr.size();
122 double dy=(py_arr[py_arr.size()-1]-py_arr[0])/py_arr.size();
123 double dz=(pz_arr[pz_arr.size()-1]-pz_arr[0])/pz_arr.size();
124 norm*=dx*dy*dz;
125
126 // Print norm
127 printf("The norm of the EMD on the cube is %e.\n",norm);
128 }
129