MADNESS 0.10.1
esolver.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30
31 $Id$
32*/
33/*
34 * File: esolver.h
35 * Author: wsttiger
36 *
37 * Created on April 17, 2009, 12:12 PM
38 */
39
40#ifndef _ESOLVER_H
41#define _ESOLVER_H
42
43typedef std::shared_ptr< WorldDCPmapInterface< Key<3> > > pmapT;
44typedef Vector<double,3> coordT;
45typedef std::shared_ptr< FunctionFunctorInterface<std::complex<double>,3> > functorT;
46typedef std::shared_ptr< FunctionFunctorInterface<double,3> > rfunctorT;
47typedef Function<std::complex<double>,3> functionT;
48typedef Function<std::complex<double>,3> cfunctionT;
49typedef Function<double,3> rfunctionT;
50typedef std::vector<functionT> vecfuncT;
51typedef std::vector<rfunctionT> rvecfuncT;
52typedef std::vector<cfunctionT> cvecfuncT;
53typedef Tensor< std::complex<double> > ctensorT;
54typedef Tensor<double> rtensorT;
55typedef FunctionFactory<std::complex<double>,3> factoryT;
56typedef FunctionFactory<double,3> rfactoryT;
57typedef SeparatedConvolution<double,3> operatorT;
58typedef std::shared_ptr<operatorT> poperatorT;
59
60void print_cube(World& world, const Function<double,3>& f, int npts)
61{
62 f.reconstruct();
63 if (world.rank() == 0) printf("\n");
64 Tensor<double> csize = FunctionDefaults<3>::get_cell_width();
65
66 for (int i = 0; i < npts; i++)
67 {
68 for (int j = 0; j < npts; j++)
69 {
70 for (int k = 0; k < npts; k++)
71 {
72 double x = (i+0.5) * (csize[0]/npts) - csize[0]/2;
73 double y = (j+0.5) * (csize[1]/npts) - csize[1]/2;
74 double z = (k+0.5) * (csize[2]/npts) - csize[2]/2;
75 coordT p(0.0);
76 p[0] = x; p[1] = y; p[2] = z;
77 if (world.rank() == 0)
78 printf("%10.2f%10.2f%10.2f%15.8f\n", x, y, z, f(p));
79 }
80 }
81 }
82}
83
84void print_cube(World& world, const Function<double,3>& f1, const Function<double,3>& f2, int npts)
85{
86 f1.reconstruct();
87 f2.reconstruct();
88 if (world.rank() == 0) printf("\n");
89 Tensor<double> csize = FunctionDefaults<3>::get_cell_width();
90
91 for (int i = 0; i < npts; i++)
92 {
93 for (int j = 0; j < npts; j++)
94 {
95 for (int k = 0; k < npts; k++)
96 {
97 double x = (i+0.5) * (csize[0]/npts) - csize[0]/2;
98 double y = (j+0.5) * (csize[1]/npts) - csize[1]/2;
99 double z = (k+0.5) * (csize[2]/npts) - csize[2]/2;
100 coordT p(0.0);
101 p[0] = x; p[1] = y; p[2] = z;
102 if (world.rank() == 0)
103 printf("%10.2f%10.2f%10.2f%15.8f%15.8f\n", x, y, z, f1(p), f2(p));
104 }
105 }
106 }
107}
108
109void print_cube(World& world, const Function<double,3>& f1, const Function<double,3>& f2,
110 const Function<double,3>& f3, int npts)
111{
112 f1.reconstruct();
113 f2.reconstruct();
114 f3.reconstruct();
115 if (world.rank() == 0) printf("\n");
116 Tensor<double> csize = FunctionDefaults<3>::get_cell_width();
117
118 for (int i = 0; i < npts; i++)
119 {
120 for (int j = 0; j < npts; j++)
121 {
122 for (int k = 0; k < npts; k++)
123 {
124 double x = (i+0.5) * (csize[0]/npts) - csize[0]/2;
125 double y = (j+0.5) * (csize[1]/npts) - csize[1]/2;
126 double z = (k+0.5) * (csize[2]/npts) - csize[2]/2;
127 coordT p(0.0);
128 p[0] = x; p[1] = y; p[2] = z;
129 if (world.rank() == 0)
130 printf("%10.2f%10.2f%10.2f%15.8f%15.8f%15.8f\n", x, y, z, f1(p), f2(p), f3(p));
131 }
132 }
133 }
134}
135
136struct KPoint
137{
138 coordT k;
139 double weight;
140 // the first wavefunction for this k-point has index begin
141 // the last wavefunction for this k-point has index end-1
142 unsigned int begin;
143 unsigned int end;
144
145 KPoint()
146 {
147 k[0] = 0.0; k[1] = 0.0; k[2] = 0.0;
148 weight = 0.0;
149 begin = -1;
150 end = -1;
151 }
152
153 KPoint(const coordT& k, const double& weight, const int& begin,
154 const int& end)
155 : k(k), weight(weight), begin(begin), end(end) {}
156
157 KPoint(const coordT& k, const double& weight)
158 : k(k), weight(weight), begin(-1), end(-1) {}
159
160 KPoint(const double& k0, const double& k1, const double& k2, const double& weight)
161 : weight(weight), begin(-1), end(-1)
162 {
163 k[0] = k0; k[1] = k1; k[2] = k2;
164 }
165
166 bool is_orb_in_kpoint(unsigned int idx)
167 {
168 return ((idx >= begin) && (idx < end)) ? true : false;
169 }
170
171 template <typename Archive>
172 void serialize(Archive& ar) {
173 ar & k & weight & begin & end;
174 }
175
176};
177
178std::istream& operator >> (std::istream& is, KPoint& kpt)
179{
180 for (unsigned int i = 0; i < kpt.k.size(); i++)
181 is >> kpt.k[i];
182 is >> kpt.weight;
183 is >> kpt.begin;
184 is >> kpt.end;
185
186 return is;
187}
188
189 //***************************************************************************
190 bool is_equal(double val1, double val2, double eps)
191 {
192 double d = fabs(val1-val2);
193 return (fabs(d) <= eps) ? true : false;
194 }
195 //***************************************************************************
196
197// //***************************************************************************
198// template <typename Q, int NDIM>
199// Function<Q,NDIM> pdiff(const Function<Q,NDIM>& f, int axis, bool fence = true)
200// {
201// Function<Q,NDIM>& g = const_cast< Function<Q,NDIM>& >(f);
202// // Check for periodic boundary conditions
203// Tensor<int> oldbc = g.get_bc();
204// Tensor<int> bc(NDIM,2);
205// bc(___) = 1;
206// g.set_bc(bc);
207// // Do calculation
208// Function<Q,NDIM> rf = diff(g,axis,fence);
209// // Restore previous boundary conditions
210// g.set_bc(oldbc);
211// return rf;
212// }
213// //***************************************************************************
214
215 //***************************************************************************
216 template <typename Q, int NDIM>
217 ctensorT kinetic_energy_matrix_slow(World& world,
218 const std::vector< Function<std::complex<Q>,NDIM> >& v,
219 const bool periodic,
220 const KPoint k = KPoint(coordT(0.0), 0.0))
221 {
222 reconstruct(world, v);
223 int n = v.size();
224 ctensorT c(n, n);
225 const std::complex<double> I = std::complex<double>(0.0, 1.0);
226 double k0 = k.k[0];
227 double k1 = k.k[1];
228 double k2 = k.k[2];
229 double ksquared = k0*k0 + k1*k1 + k2*k2;
230 if (periodic)
231 {
232 complex_derivative_3d Dx(world,0);
233 complex_derivative_3d Dy(world,1);
234 complex_derivative_3d Dz(world,2);
235 for (int i = 0; i < n; i++)
236 {
237 for (int j = 0; j <= i; j++)
238 {
239// functionT dv2_j = pdiff(pdiff(v[j], 0), 0) +
240// pdiff(pdiff(v[j], 1), 1) +
241// pdiff(pdiff(v[j], 2), 2);
242// functionT dv_j = std::complex<Q>(0.0, 2.0*k0) * pdiff(v[j], 0) +
243// std::complex<Q>(0.0, 2.0*k1) * pdiff(v[j], 1) +
244// std::complex<Q>(0.0, 2.0*k2) * pdiff(v[j], 2);
245 functionT dv2_j = Dx(Dx(v[j])) + Dy(Dy(v[j])) + Dz(Dz(v[j]));
246 functionT dv_j = std::complex<Q>(0.0, 2.0*k0) * Dx(v[j]) +
247 std::complex<Q>(0.0, 2.0*k1) * Dy(v[j]) +
248 std::complex<Q>(0.0, 2.0*k2) * Dz(v[j]);
249 functionT tmp = (ksquared*v[j]) - dv_j - dv2_j;
250 c(i, j) = inner(v[i], tmp);
251 c(j, i) = conj(c(i, j));
252 }
253 }
254 }
255 else
256 {
257 std::vector< std::shared_ptr < Derivative< std::complex<Q>,NDIM> > > gradop =
258 gradient_operator<std::complex<Q>,NDIM>(world);
259 for (int axis = 0; axis < 3; axis++)
260 {
261 vecfuncT dv = apply(world, *(gradop[axis]), v);
262 c += matrix_inner(world, dv, dv, true);
263 dv.clear(); // Allow function memory to be freed
264 }
265 }
266 return c.scale(0.5);
267 }
268
269
270 //***************************************************************************
271 template <typename Q, int NDIM>
272 ctensorT kinetic_energy_matrix(World& world,
273 const std::vector< Function<std::complex<Q>,NDIM> >& v,
274 const bool periodic,
275 const KPoint k = KPoint(coordT(0.0), 0.0))
276 {
277 const std::complex<double> I = std::complex<double>(0.0, 1.0);
278 const std::complex<double> one = std::complex<double>(1.0, 0.0);
279
280 int n = v.size();
281 ctensorT c(n, n);
282 std::vector< std::shared_ptr < Derivative< std::complex<Q>,NDIM> > > gradop =
283 gradient_operator<std::complex<Q>,NDIM>(world);
284 for (int axis = 0; axis < 3; axis++) {
285 reconstruct(world, v);
286 vecfuncT dv = apply(world, *(gradop[axis]), v);
287 if (periodic) {
288 compress(world,v);
289 compress(world,dv);
290 for (int i=0; i<n; i++) dv[i].gaxpy(one, v[i], I*k.k[axis], false);
291 world.gop.fence();
292 }
293 c += matrix_inner(world, dv, dv, true);
294 dv.clear(); // Allow function memory to be freed
295 }
296 return c.scale(0.5);
297 }
298 //***************************************************************************
299
300 //***************************************************************************
301 template <typename Q, int NDIM>
302 ctensorT kinetic_energy_matrix2(World& world,
303 const std::vector< Function<std::complex<Q>,NDIM> >& v)
304 {
305 const std::complex<double> I = std::complex<double>(0.0, 1.0);
306 const std::complex<double> one = std::complex<double>(1.0, 0.0);
307
308 int n = v.size();
309 ctensorT c(n, n);
310 std::vector< std::shared_ptr < Derivative< std::complex<Q>,NDIM> > > gradop =
311 gradient_operator<std::complex<Q>,NDIM>(world);
312 for (int axis = 0; axis < 3; axis++) {
313 reconstruct(world, v);
314 vecfuncT dv = apply(world, *(gradop[axis]), v);
315 c += matrix_inner(world, dv, dv, true);
316 dv.clear(); // Allow function memory to be freed
317 }
318 return c.scale(0.5);
319 }
320 //***************************************************************************
321
322// //***************************************************************************
323// template <typename Q, int NDIM>
324// ctensorT kinetic_energy_matrix(World& world,
325// const std::vector< Function<Q,NDIM> >& v,
326// const bool periodic,
327// const KPoint k = KPoint(coordT(0.0), 0.0))
328// {
329// reconstruct(world, v);
330// int n = v.size();
331// ctensorT c(n, n);
332// const std::complex<double> I = std::complex<double>(0.0, 1.0);
333// double k0 = k.k[0];
334// double k1 = k.k[1];
335// double k2 = k.k[2];
336// if (periodic)
337// {
338// for (int i = 0; i < n; i++)
339// {
340// functionT dv_i_0 = function_real2complex(pdiff(v[i], 0)) + I * k0 * v[i];
341// functionT dv_i_1 = function_real2complex(pdiff(v[i], 1)) + I * k1 * v[i];
342// functionT dv_i_2 = function_real2complex(pdiff(v[i], 2)) + I * k2 * v[i];
343// for (int j = 0; j <= i; j++)
344// {
345// functionT dv_j_0 = function_real2complex(pdiff(v[j], 0)) - I * k0 * v[j];
346// functionT dv_j_1 = function_real2complex(pdiff(v[j], 1)) - I * k1 * v[j];
347// functionT dv_j_2 = function_real2complex(pdiff(v[j], 2)) - I * k2 * v[j];
348// c(i, j) = inner(dv_i_0, dv_j_0) + inner(dv_i_1, dv_j_1) + inner(dv_i_2, dv_j_2);
349// c(j, i) = conj(c(i, j));
350// }
351// }
352// }
353// else
354// {
355// rtensorT r(n, n);
356// for (int axis = 0; axis < 3; axis++)
357// {
358// rvecfuncT dv = diff(world, v, axis);
359// r += matrix_inner(world, dv, dv, true);
360// dv.clear(); // Allow function memory to be freed
361// }
362// c = ctensorT(r);
363// }
364// return c.scale(0.5);
365// }
366// //***************************************************************************
367
368
369
370
371#endif /* _ESOLVER_H */
372
f
double(* f)(const coord_3d &)
Definition derivatives.cc:54
f1
double(* f1)(const coord_3d &)
Definition derivatives.cc:55
f3
double(* f3)(const coord_3d &)
Definition derivatives.cc:57
p
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
f2
double(* f2)(const coord_3d &)
Definition derivatives.cc:56
pmapT
std::shared_ptr< WorldDCPmapInterface< Key< 3 > > > pmapT
Definition esolver.h:43
cfunctionT
Function< std::complex< double >, 3 > cfunctionT
Definition esolver.h:48
operator>>
std::istream & operator>>(std::istream &is, KPoint &kpt)
Definition esolver.h:178
is_equal
bool is_equal(double val1, double val2, double eps)
Definition esolver.h:190
factoryT
FunctionFactory< std::complex< double >, 3 > factoryT
Definition esolver.h:55
print_cube
void print_cube(World &world, const Function< double, 3 > &f, int npts)
Definition esolver.h:60
ctensorT
Tensor< std::complex< double > > ctensorT
Definition esolver.h:53
coordT
Vector< double, 3 > coordT
Definition esolver.h:44
rfunctorT
std::shared_ptr< FunctionFunctorInterface< double, 3 > > rfunctorT
Definition esolver.h:46
rtensorT
Tensor< double > rtensorT
Definition esolver.h:54
rvecfuncT
std::vector< rfunctionT > rvecfuncT
Definition esolver.h:51
functionT
Function< std::complex< double >, 3 > functionT
Definition esolver.h:47
functorT
std::shared_ptr< FunctionFunctorInterface< std::complex< double >, 3 > > functorT
Definition esolver.h:45
operatorT
SeparatedConvolution< double, 3 > operatorT
Definition esolver.h:57
poperatorT
std::shared_ptr< operatorT> poperatorT
Definition esolver.h:58
rfunctionT
Function< double, 3 > rfunctionT
Definition esolver.h:49
kinetic_energy_matrix_slow
ctensorT kinetic_energy_matrix_slow(World &world, const std::vector< Function< std::complex< Q >, NDIM > > &v, const bool periodic, const KPoint k=KPoint(coordT(0.0), 0.0))
Definition esolver.h:217
cvecfuncT
std::vector< cfunctionT > cvecfuncT
Definition esolver.h:52
vecfuncT
std::vector< functionT > vecfuncT
Definition esolver.h:50
kinetic_energy_matrix2
ctensorT kinetic_energy_matrix2(World &world, const std::vector< Function< std::complex< Q >, NDIM > > &v)
Definition esolver.h:302
kinetic_energy_matrix
ctensorT kinetic_energy_matrix(World &world, const std::vector< Function< std::complex< Q >, NDIM > > &v, const bool periodic, const KPoint k=KPoint(coordT(0.0), 0.0))
Definition esolver.h:272
rfactoryT
FunctionFactory< double, 3 > rfactoryT
Definition esolver.h:56
matrix_inner
Tensor< std::complex< double > > matrix_inner(World &world, std::vector< Fcwf > &a, std::vector< Fcwf > &b)
Definition fcwf.cc:431
apply
Fcwf apply(World &world, real_convolution_3d &op, const Fcwf &psi)
Definition fcwf.cc:281
inner
std::complex< double > inner(const Fcwf &psi, const Fcwf &phi)
Definition fcwf.cc:275
gradop
std::vector< std::shared_ptr< real_derivative_3d > > gradop
Definition h2dft.cc:51
v
static const double v
Definition hatom_sf_dirac.cc:20
vecfuncT
vector< functionT > vecfuncT
Definition mcpfit.cc:51
coordT
Vector< double, 3 > coordT
Definition mcpfit.cc:48
d
static const double d
Definition nonlinschro.cc:121
c
static const double c
Definition relops.cc:10
k
static const long k
Definition rk.cc:44
KPoint
Definition esolver.h:137
KPoint::serialize
void serialize(Archive &ar)
Definition esolver.h:172
KPoint::end
unsigned int end
Definition esolver.h:143
KPoint::k
coordT k
Definition esolver.h:138
KPoint::begin
unsigned int begin
Definition esolver.h:142
KPoint::KPoint
KPoint(const coordT &k, const double &weight, const int &begin, const int &end)
Definition esolver.h:153
KPoint::KPoint
KPoint()
Definition esolver.h:145
KPoint::weight
double weight
Definition esolver.h:139
KPoint::KPoint
KPoint(const coordT &k, const double &weight)
Definition esolver.h:157
KPoint::KPoint
KPoint(const double &k0, const double &k1, const double &k2, const double &weight)
Definition esolver.h:160
KPoint::is_orb_in_kpoint
bool is_orb_in_kpoint(unsigned int idx)
Definition esolver.h:166
I
static const double_complex I
Definition tdse1d.cc:164
conj
template Tensor< double_complex > conj(const Tensor< double_complex > &t)
NDIM
static const std::size_t NDIM
Definition testpdiff.cc:42
axis
std::size_t axis
Definition testpdiff.cc:59
k0
double k0
Definition testperiodic.cc:66
k2
double k2
Definition testperiodic.cc:68
k1
double k1
Definition testperiodic.cc:67
Generated on Thu Dec 19 2024 21:22:52 for MADNESS by doxygen 1.9.8