1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
9 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
10 //
11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //////////////////////////////////////////////////////////////////////////////////////
13
14
15 #ifndef QMCPLUSPLUS_RADIALGRIDFUNCTOR_GAUSSIANBASISSET_H
16 #define QMCPLUSPLUS_RADIALGRIDFUNCTOR_GAUSSIANBASISSET_H
17 #include "hdf/hdf_archive.h"
18 #include "OhmmsData/AttributeSet.h"
19 #include <cmath>
20 #include "Message/CommOperators.h"
21 namespace qmcplusplus
22 {
23 template<class T>
24 struct GaussianCombo
25 {
26 static_assert(std::is_floating_point<T>::value, "T must be a float point type");
27 // Caution: most other code assumes value_type can only be real
28 // but maybe it can be different precision
29 // Possibly one of these types is the full precision and the other reduced precision
30 typedef T real_type;
31 real_type Y, dY, d2Y, d3Y;
32
33 struct BasicGaussian
34 {
35 real_type Sigma;
36 real_type Coeff;
37
38 real_type MinusSigma;
39 real_type CoeffP;
40 real_type CoeffPP;
41 real_type CoeffPPP1;
42 real_type CoeffPPP2;
BasicGaussianGaussianCombo::BasicGaussian43 BasicGaussian() : Sigma(1.0), Coeff(1.0) {}
44
BasicGaussianGaussianCombo::BasicGaussian45 inline BasicGaussian(real_type sig, real_type c) { reset(sig, c); }
46
resetGaussianCombo::BasicGaussian47 inline void reset(real_type sig, real_type c)
48 {
49 Sigma = sig;
50 MinusSigma = -sig;
51 Coeff = c;
52 CoeffP = -2.0 * Sigma * Coeff;
53 CoeffPP = 4.0 * Sigma * Sigma * Coeff;
54 CoeffPPP1 = 12.0 * Sigma * Sigma * Coeff;
55 CoeffPPP2 = -8.0 * Sigma * Sigma * Sigma * Coeff;
56 }
57
resetGaussianCombo::BasicGaussian58 inline void reset()
59 {
60 MinusSigma = -Sigma;
61 CoeffP = -2.0 * Sigma * Coeff;
62 CoeffPP = 4.0 * Sigma * Sigma * Coeff;
63 CoeffPPP1 = 12.0 * Sigma * Sigma * Coeff;
64 CoeffPPP2 = -8.0 * Sigma * Sigma * Sigma * Coeff;
65 }
66
setgridGaussianCombo::BasicGaussian67 inline void setgrid(real_type r) {}
68
fGaussianCombo::BasicGaussian69 inline real_type f(real_type rr) const { return Coeff * std::exp(MinusSigma * rr); }
dfGaussianCombo::BasicGaussian70 inline real_type df(real_type r, real_type rr) const { return CoeffP * r * std::exp(MinusSigma * rr); }
evaluateGaussianCombo::BasicGaussian71 inline real_type evaluate(real_type r, real_type rr, real_type& du, real_type& d2u)
72 {
73 real_type v = std::exp(MinusSigma * rr);
74 du += CoeffP * r * v;
75 d2u += (CoeffP + CoeffPP * rr) * v;
76 return Coeff * v;
77 }
evaluateGaussianCombo::BasicGaussian78 inline real_type evaluate(real_type r, real_type rr, real_type& du, real_type& d2u, real_type& d3u)
79 {
80 real_type v = std::exp(MinusSigma * rr);
81 du += CoeffP * r * v;
82 d2u += (CoeffP + CoeffPP * rr) * v;
83 d3u += (CoeffPPP1 * r + CoeffPPP2 * r * rr) * v;
84 return Coeff * v;
85 }
86 };
87
88 ///Boolean
89 bool Normalized;
90 real_type L;
91 real_type NormL;
92 real_type NormPow;
93 std::string nodeName;
94 std::string expName;
95 std::string coeffName;
96 std::vector<BasicGaussian> gset;
97
98 explicit GaussianCombo(int l = 0,
99 bool normalized = false,
100 const char* node_name = "radfunc",
101 const char* exp_name = "exponent",
102 const char* c_name = "contraction");
103
104 void reset();
105
106 /** return the number Gaussians
107 */
sizeGaussianCombo108 inline int size() const { return gset.size(); }
109
fGaussianCombo110 inline real_type f(real_type r)
111 {
112 real_type res = 0;
113 real_type r2 = r * r;
114 typename std::vector<BasicGaussian>::const_iterator it(gset.begin());
115 typename std::vector<BasicGaussian>::const_iterator it_end(gset.end());
116 while (it != it_end)
117 {
118 res += (*it).f(r2);
119 ++it;
120 }
121 return res;
122 }
dfGaussianCombo123 inline real_type df(real_type r)
124 {
125 real_type res = 0;
126 real_type r2 = r * r;
127 typename std::vector<BasicGaussian>::const_iterator it(gset.begin());
128 typename std::vector<BasicGaussian>::const_iterator it_end(gset.end());
129 while (it != it_end)
130 {
131 res += (*it).df(r, r2);
132 ++it;
133 }
134 return res;
135 }
136
evaluateGaussianCombo137 inline real_type evaluate(real_type r, real_type rinv)
138 {
139 Y = 0.0;
140 real_type rr = r * r;
141 typename std::vector<BasicGaussian>::iterator it(gset.begin()), it_end(gset.end());
142 while (it != it_end)
143 {
144 Y += (*it).f(rr);
145 ++it;
146 }
147 return Y;
148 }
149
evaluateAllGaussianCombo150 inline void evaluateAll(real_type r, real_type rinv)
151 {
152 Y = 0.0;
153 dY = 0.0;
154 d2Y = 0.0;
155 real_type rr = r * r;
156 typename std::vector<BasicGaussian>::iterator it(gset.begin()), it_end(gset.end());
157 while (it != it_end)
158 {
159 Y += (*it).evaluate(r, rr, dY, d2Y);
160 ++it;
161 }
162 }
163
evaluateWithThirdDerivGaussianCombo164 inline void evaluateWithThirdDeriv(real_type r, real_type rinv)
165 {
166 Y = 0.0;
167 dY = 0.0;
168 d2Y = 0.0, d3Y = 0.0;
169 real_type rr = r * r;
170 typename std::vector<BasicGaussian>::iterator it(gset.begin()), it_end(gset.end());
171 while (it != it_end)
172 {
173 Y += (*it).evaluate(r, rr, dY, d2Y, d3Y);
174 ++it;
175 }
176 }
177
178 //inline real_type evaluate(real_type r, real_type rinv, real_type& drnl, real_type& d2rnl) {
179 // Y=0.0;drnl=0.0;d2rnl=0.0;
180 // real_type du, d2u;
181 // typename std::vector<BasicGaussian>::iterator
182 // it(sset.begin()),it_end(sset.end());
183 // while(it != it_end) {
184 // Y+=(*it).evaluate(r,rinv,du,d2u); dY+=du; d2Y+=d2u;
185 // ++it;
186 // }
187 // return Y;
188 //}
189
190
191 bool put(xmlNodePtr cur);
192
193 void addGaussian(real_type c, real_type alpha);
194
195 bool putBasisGroup(xmlNodePtr cur);
196 bool putBasisGroupH5(hdf_archive& hin);
197
198 /** double factorial of num
199 * @param num integer to be factored
200 * @return num!!
201 *
202 * \if num == odd,
203 * \f$ num!! = 1\cdot 3\cdot ... \cdot num-2 \cdot num\f$
204 * \else num == even,
205 * \f$ num!! = 2\cdot 4\cdot ... \cdot num-2 \cdot num\f$
206 */
DFactorialGaussianCombo207 int DFactorial(int num) { return (num < 2) ? 1 : num * DFactorial(num - 2); }
208 };
209
210 template<class T>
GaussianCombo(int l,bool normalized,const char * node_name,const char * exp_name,const char * c_name)211 GaussianCombo<T>::GaussianCombo(int l, bool normalized, const char* node_name, const char* exp_name, const char* c_name)
212 : Normalized(normalized), nodeName(node_name), expName(exp_name), coeffName(c_name)
213 {
214 L = static_cast<real_type>(l);
215 //Everything related to L goes to NormL and NormPow
216 const real_type pi = 4.0 * std::atan(1.0);
217 NormL =
218 std::pow(2, L + 1) * std::sqrt(2.0 / static_cast<real_type>(DFactorial(2 * l + 1))) * std::pow(2.0 / pi, 0.25);
219 NormPow = 0.5 * (L + 1.0) + 0.25;
220 }
221
222 template<class T>
put(xmlNodePtr cur)223 bool GaussianCombo<T>::put(xmlNodePtr cur)
224 {
225 real_type alpha(1.0), c(1.0);
226 OhmmsAttributeSet radAttrib;
227 radAttrib.add(alpha, expName);
228 radAttrib.add(c, coeffName);
229 radAttrib.put(cur);
230 addGaussian(c, alpha);
231 return true;
232 }
233
234 template<class T>
addGaussian(real_type c,real_type alpha)235 void GaussianCombo<T>::addGaussian(real_type c, real_type alpha)
236 {
237 if (!Normalized)
238 {
239 c *= NormL * std::pow(alpha, NormPow);
240 }
241 // LOGMSG(" Gaussian exponent = " << alpha
242 // << "\n contraction=" << c0 << " normalized contraction = " << c)
243 gset.push_back(BasicGaussian(alpha, c));
244 }
245
246 template<class T>
reset()247 void GaussianCombo<T>::reset()
248 {
249 for (int i = 0; i < gset.size(); i++)
250 gset[i].reset();
251 }
252
253 template<class T>
putBasisGroup(xmlNodePtr cur)254 bool GaussianCombo<T>::putBasisGroup(xmlNodePtr cur)
255 {
256 cur = cur->children;
257 while (cur != NULL)
258 {
259 std::string cname((const char*)cur->name);
260 if (cname == "radfunc")
261 {
262 put(cur);
263 }
264 cur = cur->next;
265 }
266 reset();
267 return true;
268 }
269
270 template<class T>
putBasisGroupH5(hdf_archive & hin)271 bool GaussianCombo<T>::putBasisGroupH5(hdf_archive& hin)
272 {
273 int NbRadFunc(0);
274 if (hin.myComm->rank() == 0)
275 {
276 hin.read(NbRadFunc, "NbRadFunc");
277 hin.push("radfunctions");
278 }
279 hin.myComm->bcast(NbRadFunc);
280
281 for (int i = 0; i < NbRadFunc; i++)
282 {
283 real_type alpha(1.0), c(1.0);
284 std::stringstream tempdata;
285 std::string dataradID0 = "DataRad", dataradID;
286 tempdata << dataradID0 << i;
287 dataradID = tempdata.str();
288
289 if (hin.myComm->rank() == 0)
290 {
291 hin.push(dataradID.c_str());
292 hin.read(alpha, "exponent");
293 hin.read(c, "contraction");
294 }
295
296 hin.myComm->bcast(alpha);
297 hin.myComm->bcast(c);
298
299 if (!Normalized)
300 c *= NormL * std::pow(alpha, NormPow);
301 // LOGMSG(" Gaussian exponent = " << alpha
302 // << "\n contraction=" << c0 << " nomralized contraction = " << c)
303 gset.push_back(BasicGaussian(alpha, c));
304 if (hin.myComm->rank() == 0)
305 hin.pop();
306 }
307 reset();
308 if (hin.myComm->rank() == 0)
309 hin.pop();
310
311 return true;
312 }
313 } // namespace qmcplusplus
314 #endif
315