MADNESS  0.10.1
electronicstructureapp.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  * electronicstructureapp.h
35  *
36  * Created on: Nov 5, 2008
37  * Author: eh7
38  */
39 
40 #ifndef ELECTRONICSTRUCTUREAPP_H_
41 #define ELECTRONICSTRUCTUREAPP_H_
42 
43 #include <madness/mra/mra.h>
44 #include <madness/misc/ran.h>
46 //#include "poperator.h"
47 #include "libxc.h"
48 #include "complexfun.h"
49 #include "esolver.h"
50 
51 class LevelPmap : public WorldDCPmapInterface< Key<3> > {
52 private:
53  const int nproc;
54 public:
55  LevelPmap() : nproc(0) {};
56 
57  LevelPmap(World& world) : nproc(world.nproc()) {}
58 
59  /// Find the owner of a given key
60  ProcessID owner(const Key<3>& key) const {
61  Level n = key.level();
62  if (n == 0) return 0;
63  hashT hash;
64  if (n <= 3 || (n&0x1)) hash = key.hash();
65  else hash = key.parent().hash();
66  //hashT hash = key.hash();
67  return hash%nproc;
68  }
69 };
70 
71 class MolecularPotentialFunctor : public FunctionFunctorInterface<double,3> {
72 private:
74 public:
76  : _mentity(mentity)
77  {}
78 
79  double operator()(const coordT& x) const
80  {
81  return _mentity.nuclear_attraction_potential(x[0], x[1], x[2]);
82  }
83 };
84 
86 private:
88  const double R;
89  const bool periodic;
90  const std::vector<coordT> _specialpts;
91 public:
93  const bool& periodic, const std::vector<coordT>& specialpts)
94  : _mentity(mentity), R(R), periodic(periodic), _specialpts(specialpts) {
95  }
96 
97  virtual std::vector<coordT> special_points() const
98  {
99  return _specialpts;
100  }
101 
103  {
104  return 10;
105  }
106 
107  double operator()(const coordT& x) const
108  {
109  //double big = 0.5*R + 6.0*_mentity.smallest_length_scale();
110  double big = 2*R + 6.0*_mentity.smallest_length_scale();
111  // Only one contribution at any point due to the short
112  // range of the nuclear charge density
113  //printf("big: %10.8f\n\n", big);
114  double value = 0.0;
115  if (periodic)
116  {
117  for (int xr = -1; xr <= 1; xr += 1)
118  {
119  double xx = x[0] + xr*R;
120  //printf("x[0]: %10.8f xx: %10.8f\n", x[0], xx);
121  if (xx < big && xx > -big)
122  {
123  for (int yr = -1; yr <= 1; yr += 1)
124  {
125  double yy = x[1] + yr*R;
126  //printf("y[0]: %10.8f yy: %10.8f\n", x[1], yy);
127  if (yy < big && yy > -big)
128  {
129  for (int zr = -1; zr <= 1; zr += 1)
130  {
131  double zz = x[2] + zr*R;
132  //printf("z[0]: %10.8f zz: %10.8f\n", x[2], zz);
133  if (zz < big && zz > -big)
134  {
135  double t1 = _mentity.nuclear_charge_density(xx, yy, zz);
136  value += t1;
137  //printf("t1: %10.8f value: %10.8f\n", t1, value);
138  }
139  }
140  }
141  }
142  }
143  }
144  }
145  else
146  {
147  value = _mentity.nuclear_charge_density(x[0], x[1], x[2]);
148  }
149  return value;
150  }
151 };
152 
153 #define NTRANS 8
154 
155 class AtomicBasisFunctor : public FunctionFunctorInterface<std::complex<double>,3> {
156 private:
158  const double R;
159  const double rangesq;
160  const bool periodic;
161  const KPoint kpt;
163  std::vector<coordT> _specialpts;
164 // const int NTRANS = 2;
165  Vector<std::complex<double>,2*NTRANS+1> tx;
166  Vector<std::complex<double>,2*NTRANS+1> ty;
167  Vector<std::complex<double>,2*NTRANS+1> tz;
168 public:
170  bool periodic, const KPoint kpt)
172  {
173  double x, y, z;
174  aofunc.get_coords(x,y,z);
175 
176  r[0]=x; r[1]=y; r[2]=z;
177  _specialpts=std::vector<coordT>(1,r);
178 
179  for (int ir = -NTRANS; ir <= NTRANS; ir += 1)
180  {
181  tx[ir+NTRANS] = exp(std::complex<double>(0.0, kpt.k[0]*ir * R));
182  ty[ir+NTRANS] = exp(std::complex<double>(0.0, kpt.k[1]*ir * R));
183  tz[ir+NTRANS] = exp(std::complex<double>(0.0, kpt.k[2]*ir * R));
184  }
185 }
186 
187  virtual std::vector<coordT> special_points() const
188  {
189  return _specialpts;
190  }
191 
192  std::complex<double> operator()(const coordT& x) const
193  {
194  std::complex<double> value = 0.0;
195  if (periodic) {
196  for (int xx=-NTRANS; xx<=NTRANS; xx++) {
197  const double xxR = xx*R + x[0] -r[0];
198  const double xxRsq = xxR*xxR;
199  if (xxRsq < rangesq) {
200  for (int yy=-NTRANS; yy<=NTRANS; yy++) {
201  const double yyR = yy*R + x[1] - r[1];
202  const double yyRsq = xxRsq + yyR*yyR;
203  if (yyRsq < rangesq) {
204  for (int zz=-NTRANS; zz<=NTRANS; zz++) {
205  double ao = aofunc(xx*R+x[0], yy*R+x[1], zz*R+x[2]);
206  if (fabs(ao) > 1e-8) {
207  std::complex<double> t1 = tx[xx+NTRANS]*ty[yy+NTRANS]*tz[zz+NTRANS];
208  double kx0 = kpt.k[0] * x[0];
209  double kx1 = kpt.k[1] * x[1];
210  double kx2 = kpt.k[2] * x[2];
211  std::complex<double> t2 = exp(std::complex<double>(0.0, -kx0 - kx1 - kx2));
212 // value += t1 * t2 * ao;
213 // value += t1 * ao;
214  value += ao;
215  }
216  }
217  }
218  }
219  }
220  }
221  }
222  else {
223  value = aofunc(x[0], x[1], x[2]);
224  }
225  return value;
226  }
227 };
228 
229 //std::complex<double> operator()(const coordT& x) const
230 //{
231 // std::complex<double> value = 0.0;
232 // if (periodic) {
233 // for (int xx=-NTRANS; xx<=NTRANS; xx++) {
234 // const double xxR = x[0] - r[0] - xx*R;
235 // const double xxRsq = xxR*xxR;
236 // if (xxRsq < rangesq) {
237 // for (int yy=-NTRANS; yy<=NTRANS; yy++) {
238 // const double yyR = x[1] - r[1] - yy*R;
239 // const double yyRsq = xxRsq + yyR*yyR;
240 // if (yyRsq < rangesq) {
241 // for (int zz=-NTRANS; zz<=NTRANS; zz++) {
242 // const double zzR = x[2] - r[2] - zz*R;
243 // const double zzRsq = yyRsq + zzR*zzR;
244 // if (zzRsq < rangesq) {
245 // double ao = aofunc(x[0]-xx*R, x[1]-yy*R, x[2]-zz*R);
246 // if (fabs(ao) > 1e-8) {
247 // std::complex<double> t1 = tx[xx+NTRANS]*ty[yy+NTRANS]*tz[zz+NTRANS];
248 // double kx0 = kpt.k[0] * x[0];
249 // double kx1 = kpt.k[1] * x[1];
250 // double kx2 = kpt.k[2] * x[2];
251 // std::complex<double> t2 = exp(std::complex<double>(0.0, -kx0 - kx1 - kx2));
252 // value += t1 * t2 * ao;
253 // // value += t1 * ao;
254 // }
255 // }
256 // }
257 // }
258 // }
259 // }
260 // }
261 // }
262 // else {
263 // value = aofunc(x[0], x[1], x[2]);
264 // }
265 // return value;
266 //}
267 //};
268 
269 double rsquared(const coordT& r) {
270  return r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
271 }
272 
273 
274 #endif /* ELECTRONICSTRUCTUREAPP_H_ */
Used to represent one basis function from a shell on a specific center.
Definition: apps/periodic_old/molecularbasis.h:361
void get_coords(double &x, double &y, double &z) const
Definition: apps/periodic_old/molecularbasis.h:408
Definition: preal.cc:77
Vector< std::complex< double >, 2 *NTRANS+1 > tx
Definition: electronicstructureapp.h:165
Vector< std::complex< double >, 2 *NTRANS+1 > tz
Definition: electronicstructureapp.h:167
std::complex< double > operator()(const coordT &x) const
Definition: electronicstructureapp.h:192
const KPoint kpt
Definition: electronicstructureapp.h:161
coordT r
Definition: electronicstructureapp.h:162
virtual std::vector< coordT > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: electronicstructureapp.h:187
AtomicBasisFunctor(const AtomicBasisFunction &aofunc, double R, bool periodic, const KPoint kpt)
Definition: electronicstructureapp.h:169
const AtomicBasisFunction aofunc
Definition: preal.cc:79
const double R
Definition: electronicstructureapp.h:158
const double rangesq
Definition: electronicstructureapp.h:159
std::vector< coordT > _specialpts
Definition: electronicstructureapp.h:163
Vector< std::complex< double >, 2 *NTRANS+1 > ty
Definition: electronicstructureapp.h:166
const bool periodic
Definition: electronicstructureapp.h:160
Definition: electronicstructureapp.h:51
ProcessID owner(const Key< 3 > &key) const
Find the owner of a given key.
Definition: electronicstructureapp.h:60
LevelPmap()
Definition: electronicstructureapp.h:55
LevelPmap(World &world)
Definition: electronicstructureapp.h:57
const int nproc
Definition: electronicstructureapp.h:53
Definition: mentity.h:95
double nuclear_attraction_potential(double x, double y, double z) const
Definition: mentity.cc:452
double smallest_length_scale() const
Definition: mentity.cc:403
double nuclear_charge_density(double x, double y, double z) const
Definition: mentity.cc:469
Definition: electronicstructureapp.h:85
virtual Level special_level()
Override this change level refinement for special points (default is 6)
Definition: electronicstructureapp.h:102
const MolecularEntity & _mentity
Definition: electronicstructureapp.h:87
double operator()(const coordT &x) const
Definition: electronicstructureapp.h:107
const double R
Definition: electronicstructureapp.h:88
const std::vector< coordT > _specialpts
Definition: electronicstructureapp.h:90
const bool periodic
Definition: electronicstructureapp.h:89
virtual std::vector< coordT > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: electronicstructureapp.h:97
MolecularNuclearChargeDensityFunctor(const MolecularEntity &mentity, const double &R, const bool &periodic, const std::vector< coordT > &specialpts)
Definition: electronicstructureapp.h:92
Definition: preal.cc:121
MolecularPotentialFunctor(const MolecularEntity &mentity)
Definition: electronicstructureapp.h:75
const MolecularEntity & _mentity
Definition: electronicstructureapp.h:73
double operator()(const coordT &x) const
Definition: electronicstructureapp.h:79
FunctionFunctorInterface()
Definition: function_interface.h:77
double rsquared(const coordT &r)
Definition: electronicstructureapp.h:269
#define NTRANS
Definition: electronicstructureapp.h:153
Vector< double, 3 > coordT
Definition: mcpfit.cc:48
Main include file for MADNESS and defines Function interface.
int Level
Definition: key.h:55
std::size_t hashT
The hash value type.
Definition: worldhash.h:145
Definition: esolver.h:137
coordT k
Definition: esolver.h:138
Definition: test_dc.cc:47
Key parent() const
Definition: test_tree.cc:68
hashT hash() const
Definition: test_dc.cc:54
void e()
Definition: test_sig.cc:75
int ProcessID
Used to clearly identify process number/rank.
Definition: worldtypes.h:43