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) 2019 QMCPACK developers.
6 //
7 // File developed by: Ken Esler, kpesler@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 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //                    Luke Shulenburger, lshulen@sandia.gov, Sandia National Laboratories
13 //
14 // File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
15 //////////////////////////////////////////////////////////////////////////////////////
16 
17 
18 /*
19  DO NOT MAKE PERMANENT EDITS IN THIS FILE
20  This file is generated from src/QMCWavefunctions/Jastrow/codegen/user_jastrow.py and UserFunctor.h.in
21 
22  To make changes, edit UserFunctor.h.in and rerun user_jastrow.py.
23 */
24 
25 
26 /** @file UserFunctor.h
27  * @brief User-defined functor
28  */
29 #ifndef QMCPLUSPLUS_USERFUNCTOR_H
30 #define QMCPLUSPLUS_USERFUNCTOR_H
31 #include "Numerics/OptimizableFunctorBase.h"
32 #include "OhmmsData/AttributeSet.h"
33 #include <cmath>
34 // #include <vector>
35 #include "OhmmsPETE/TinyVector.h"
36 
37 
38 namespace qmcplusplus
39 {
40 /** Implements the function
41  *  f = A*r/(B*r + 1) - A/B
42 
43  *
44  */
45 template<class T>
46 struct UserFunctor : public OptimizableFunctorBase
47 {
48 
49 
50   /// Is optimizable
51   bool Opt_A;
52   /// Value
53   real_type A;
54   /// Id in XML input
55   std::string ID_A;
56 
57 
58   /// Is optimizable
59   bool Opt_B;
60   /// Value
61   real_type B;
62   /// Id in XML input
63   std::string ID_B;
64 
65 
66 
67   ///default constructor
UserFunctorUserFunctor68   UserFunctor() { reset(); }
69 
70 // void setCusp(real_type cusp)
71 
setCuspUserFunctor72   void setCusp(real_type cusp)
73   {
74     A     = cusp;
75     Opt_A = false;
76     reset();
77   }
78 
79 
makeCloneUserFunctor80   OptimizableFunctorBase* makeClone() const { return new UserFunctor(*this); }
81 
resetUserFunctor82   void reset()
83   {
84   }
85 
86 
87 
88 // inline real_type evaluate(real_type r) const
89 
evaluateUserFunctor90   inline real_type evaluate(real_type r) const {
91    return A*r/(B*r + 1) - A/B;
92   }
93 
94 
95 // const inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) const
96 
evaluateUserFunctor97   inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) const
98   {
99     dudr        = -A*B*r/((B*r + 1)*(B*r + 1)) + A/(B*r + 1);
100     d2udr2      = 2*A*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1));
101     return A*r/(B*r + 1) - A/B;
102   }
103 
104 
105 
106 // inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3) const
107 
evaluateUserFunctor108   inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3) const
109   {
110     dudr        = -A*B*r/((B*r + 1)*(B*r + 1)) + A/(B*r + 1);
111     d2udr2      = 2*A*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1));
112     d3udr3      = 6*A*((B)*(B))*(-B*r/(B*r + 1) + 1)/std::pow(B*r + 1, 3);
113     return A*r/(B*r + 1) - A/B;
114   }
115 
116 
117 
118 
evaluateVUserFunctor119   inline real_type evaluateV(const int iat,
120                              const int iStart,
121                              const int iEnd,
122                              const T* restrict _distArray,
123                              T* restrict distArrayCompressed) const
124   {
125   // specialized evaluation loop?
126     real_type sum(0);
127     for (int idx = iStart; idx < iEnd; idx++)
128       if (idx != iat)
129         sum += evaluate(_distArray[idx]);
130     return sum;
131   }
132 
evaluateVGLUserFunctor133   inline void evaluateVGL(const int iat,
134                           const int iStart,
135                           const int iEnd,
136                           const T* distArray,
137                           T* restrict valArray,
138                           T* restrict gradArray,
139                           T* restrict laplArray,
140                           T* restrict distArrayCompressed,
141                           int* restrict distIndices) const
142   {
143     // specialized evaluation loop?
144     for (int idx = iStart; idx < iEnd; idx++)
145     {
146       valArray[idx] = evaluate(distArray[idx], gradArray[idx], laplArray[idx]);
147       gradArray[idx] /= distArray[idx];
148     }
149     if (iat >= iStart && iat < iEnd)
150       valArray[iat] = gradArray[iat] = laplArray[iat] = T(0);
151   }
152 
fUserFunctor153   inline real_type f(real_type r) { return evaluate(r); }
154 
dfUserFunctor155   inline real_type df(real_type r)
156   {
157     real_type dudr, d2udr2;
158     real_type res = evaluate(r, dudr, d2udr2);
159     return dudr;
160   }
161 
162 // inline bool evaluateDerivatives(real_type r, std::vector<TinyVector<real_type, 3>>& derivs)
163 
evaluateDerivativesUserFunctor164   inline bool evaluateDerivatives(real_type r, std::vector<TinyVector<real_type, 3>>& derivs)
165   {
166     int i = 0;
167 
168     if (Opt_A)
169     {
170       derivs[i][0] = r/(B*r + 1) - 1/B;
171       derivs[i][1] = -B*r/((B*r + 1)*(B*r + 1)) + 1.0/(B*r + 1);
172       derivs[i][2] = 2*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1));
173       ++i;
174     }
175 
176     if (Opt_B)
177     {
178       derivs[i][0] = -A*((r)*(r))/((B*r + 1)*(B*r + 1)) + A/((B)*(B));
179       derivs[i][1] = 2*A*B*((r)*(r))/std::pow(B*r + 1, 3) - 2*A*r/((B*r + 1)*(B*r + 1));
180       derivs[i][2] = -4*A*B*r*(B*r/(B*r + 1) - 1)/std::pow(B*r + 1, 3) + 2*A*B*(-B*((r)*(r))/((B*r + 1)*(B*r + 1)) + r/(B*r + 1))/((B*r + 1)*(B*r + 1)) + 2*A*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1));
181       ++i;
182     }
183 
184     return true;
185   }
186 
187 
188 // inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs)
189 
190   /// compute derivatives with respect to variational parameters
evaluateDerivativesUserFunctor191   inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs)
192   {
193     int i = 0;
194 
195     if (Opt_A)
196     {
197       derivs[i] = r/(B*r + 1) - 1/B;
198       ++i;
199     }
200 
201     if (Opt_B)
202     {
203       derivs[i] = -A*((r)*(r))/((B*r + 1)*(B*r + 1)) + A/((B)*(B));
204       ++i;
205     }
206 
207 
208     return true;
209   }
210 
211 
212 //  bool put(xmlNodePtr cur)
213 
putUserFunctor214   bool put(xmlNodePtr cur)
215   {
216     cur = cur->xmlChildrenNode;
217     while (cur != NULL)
218     {
219       std::string cname((const char*)(cur->name));
220       if (cname == "var") //only accept var
221       {
222         std::string id_in;
223         std::string p_name;
224         OhmmsAttributeSet rAttrib;
225         rAttrib.add(id_in, "id");
226         rAttrib.add(p_name, "name");
227         rAttrib.put(cur);
228 
229         if (p_name == "A")
230         {
231           ID_A = id_in;
232           putContent(A, cur);
233           Opt_A = true;
234         }
235 
236         if (p_name == "B")
237         {
238           ID_B = id_in;
239           putContent(B, cur);
240           Opt_B = true;
241         }
242 
243       }
244       cur = cur->next;
245     }
246     reset();
247     myVars.clear();
248 
249     if (Opt_A)
250       myVars.insert(ID_A, A, Opt_A, optimize::OTHER_P);
251 
252     if (Opt_B)
253       myVars.insert(ID_B, B, Opt_B, optimize::OTHER_P);
254 
255     return true;
256   }
257 
258 
checkInVariablesUserFunctor259   void checkInVariables(opt_variables_type& active)
260   {
261     active.insertFrom(myVars);
262     //myVars.print(std::cout);
263   }
264 
checkOutVariablesUserFunctor265   void checkOutVariables(const opt_variables_type& active)
266   {
267     myVars.getIndex(active);
268     //myVars.print(std::cout);
269   }
270 
271 //void resetParameters(const opt_variables_type& active)
272 
resetParametersUserFunctor273   void resetParameters(const opt_variables_type& active)
274   {
275     if (myVars.size())
276     {
277       int ia = myVars.where(0);
278       if (ia > -1)
279       {
280         int i = 0;
281 
282         if (Opt_A)
283           A = std::real(myVars[i++] = active[ia++]);
284 
285         if (Opt_B)
286           B = std::real(myVars[i++] = active[ia++]);
287 
288       }
289       reset();
290     }
291   }
292 
293 };
294 
295 
296 
297 
298 } // namespace qmcplusplus
299 #endif
300