MADNESS  version 0.9
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 
43 typedef std::shared_ptr< WorldDCPmapInterface< Key<3> > > pmapT;
44 typedef Vector<double,3> coordT;
45 typedef std::shared_ptr< FunctionFunctorInterface<std::complex<double>,3> > functorT;
46 typedef std::shared_ptr< FunctionFunctorInterface<double,3> > rfunctorT;
47 typedef Function<std::complex<double>,3> functionT;
48 typedef Function<std::complex<double>,3> cfunctionT;
49 typedef Function<double,3> rfunctionT;
50 typedef std::vector<functionT> vecfuncT;
51 typedef std::vector<rfunctionT> rvecfuncT;
52 typedef std::vector<cfunctionT> cvecfuncT;
53 typedef Tensor< std::complex<double> > ctensorT;
54 typedef Tensor<double> rtensorT;
55 typedef FunctionFactory<std::complex<double>,3> factoryT;
56 typedef FunctionFactory<double,3> rfactoryT;
57 typedef SeparatedConvolution<double,3> operatorT;
58 typedef std::shared_ptr<operatorT> poperatorT;
59 
60 void 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 
84 void 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 
109 void 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 
136 struct 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 
178 std::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 
coordT
Vector< double, 3 > coordT
Definition: esolver.h:44
pmapT
std::shared_ptr< WorldDCPmapInterface< Key< 3 > > > pmapT
Definition: esolver.h:43
madness::tr1::shptr::shared_ptr
Definition: shared_ptr_bits.h:38
madness::reconstruct
void reconstruct(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Reconstruct a vector of functions.
Definition: vmra.h:149
std::tr1::f3
const T1 const T2 const T3 & f3
Definition: gtest-tuple.h:686
madness::gaxpy
void gaxpy(World &world, Q alpha, std::vector< Function< T, NDIM > > &a, Q beta, const std::vector< Function< R, NDIM > > &b, bool fence=true)
Generalized A*X+Y for vectors of functions -— a[i] = alpha*a[i] + beta*b[i].
Definition: vmra.h:680
rvecfuncT
std::vector< rfunctionT > rvecfuncT
Definition: esolver.h:51
print_cube
void print_cube(World &world, const Function< double, 3 > &f, int npts)
Definition: esolver.h:60
KPoint::begin
unsigned int begin
Definition: esolver.h:142
cfunctionT
Function< std::complex< double >, 3 > cfunctionT
Definition: esolver.h:48
NDIM
const int NDIM
Definition: tdse1.cc:44
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
kinetic_energy_matrix2
ctensorT kinetic_energy_matrix2(World &world, const std::vector< Function< std::complex< Q >, NDIM > > &v)
Definition: esolver.h:302
vecfuncT
std::vector< functionT > vecfuncT
Definition: esolver.h:50
KPoint::weight
double weight
Definition: esolver.h:139
rfunctionT
Function< double, 3 > rfunctionT
Definition: esolver.h:49
rfactoryT
FunctionFactory< double, 3 > rfactoryT
Definition: esolver.h:56
madness::f
NDIM & f
Definition: mra.h:2179
functionT
Function< double, 3 > functionT
Definition: DFcode/corepotential.cc:57
cvecfuncT
std::vector< cfunctionT > cvecfuncT
Definition: esolver.h:52
KPoint
Definition: esolver.h:136
KPoint::serialize
void serialize(Archive &ar)
Definition: esolver.h:172
madness::inner
T inner(const vecfunc< T, NDIM > &a, const vecfunc< T, NDIM > &b)
the non-linear solver requires an inner product
Definition: nemo.h:112
factoryT
FunctionFactory< std::complex< double >, 3 > factoryT
Definition: esolver.h:55
KPoint::KPoint
KPoint()
Definition: esolver.h:145
k
const int k
Definition: dielectric.cc:184
one
const complexd one(1, 0)
coordT
Vector< double, 3 > coordT
Definition: DFcode/corepotential.cc:55
gradop
std::vector< std::shared_ptr< real_derivative_3d > > gradop
Definition: h2dft.cc:52
KPoint::k
coordT k
Definition: esolver.h:138
KPoint::end
unsigned int end
Definition: esolver.h:143
operator>>
std::istream & operator>>(std::istream &is, KPoint &kpt)
Definition: esolver.h:178
rfunctorT
std::shared_ptr< FunctionFunctorInterface< double, 3 > > rfunctorT
Definition: esolver.h:46
KPoint::KPoint
KPoint(const coordT &k, const double &weight)
Definition: esolver.h:157
I
const complexd I(0, 1)
ctensorT
Tensor< std::complex< double > > ctensorT
Definition: esolver.h:53
poperatorT
std::shared_ptr< operatorT > poperatorT
Definition: esolver.h:58
functionT
Function< std::complex< double >, 3 > functionT
Definition: esolver.h:47
rtensorT
Tensor< double > rtensorT
Definition: esolver.h:54
KPoint::KPoint
KPoint(const double &k0, const double &k1, const double &k2, const double &weight)
Definition: esolver.h:160
std::tr1::f1
const T1 & f1
Definition: gtest-tuple.h:680
madness::complex_derivative_3d
Derivative< double_complex, 3 > complex_derivative_3d
Definition: functypedefs.h:177
std::tr1::f2
const T1 const T2 & f2
Definition: gtest-tuple.h:680
is_equal
bool is_equal(double val1, double val2, double eps)
Definition: esolver.h:190
mpfr::fabs
const mpreal fabs(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2187
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
KPoint::KPoint
KPoint(const coordT &k, const double &weight, const int &begin, const int &end)
Definition: esolver.h:153
madness::apply
Function< T, NDIM > apply(const Derivative< T, NDIM > &D, const Function< T, NDIM > &f, bool fence=true)
Applies derivative operator to function (for syntactic equivalence to integral operator apply) ...
Definition: derivative.h:613
c
const double c
Definition: gfit.cc:200
operatorT
SeparatedConvolution< double, 3 > operatorT
Definition: esolver.h:57
KPoint::is_orb_in_kpoint
bool is_orb_in_kpoint(unsigned int idx)
Definition: esolver.h:166
matrix_inner
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition: DFcode/distpm.cc:38
functorT
std::shared_ptr< FunctionFunctorInterface< std::complex< double >, 3 > > functorT
Definition: esolver.h:45
madness::conj
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence...
Definition: mra.h:1879
madness::compress
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:130
vecfuncT
std::vector< Function< double, 3 > > vecfuncT
Definition: tdhf_CIS.h:75
Generated on Fri Sep 26 2014 10:38:50 for MADNESS by   doxygen 1.8.8