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
51class LevelPmap : public WorldDCPmapInterface< Key<3> > {
52private:
53 const int nproc;
54public:
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
72private:
74public:
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
86private:
88 const double R;
89 const bool periodic;
90 const std::vector<coordT> _specialpts;
91public:
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
155class AtomicBasisFunctor : public FunctionFunctorInterface<std::complex<double>,3> {
156private:
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;
168public:
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
269double 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
const KPoint kpt
Definition electronicstructureapp.h:161
virtual std::vector< coordT > special_points() const
Override this to return list of special points to be refined more deeply.
Definition electronicstructureapp.h:187
coordT r
Definition electronicstructureapp.h:162
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
std::complex< double > operator()(const coordT &x) const
Definition electronicstructureapp.h:192
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
virtual std::vector< coordT > special_points() const
Override this to return list of special points to be refined more deeply.
Definition electronicstructureapp.h:97
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
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
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
FunctionFunctorInterface()
Definition function_interface.h:77
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:66
Level level() const
Definition key.h:159
hashT hash() const
Definition key.h:148
Key parent(int generation=1) const
Returns the key of the parent.
Definition key.h:187
Interface to be provided by any process map.
Definition worlddc.h:82
A parallel world class.
Definition world.h:132
double rsquared(const coordT &r)
Definition electronicstructureapp.h:269
#define NTRANS
Definition electronicstructureapp.h:153
Vector< double, 3 > coordT
Definition ewald.cc:9
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
void e()
Definition test_sig.cc:75
int ProcessID
Used to clearly identify process number/rank.
Definition worldtypes.h:43