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 //                    Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
9 //                    Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
10 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
11 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //
13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15 
16 
17 #ifndef QMCPLUSPLUS_CUBICFUNCTORSFORJASTROW_H
18 #define QMCPLUSPLUS_CUBICFUNCTORSFORJASTROW_H
19 #include "Numerics/OneDimGridBase.h"
20 #include "Numerics/CubicBspline.h"
21 #include "Numerics/OptimizableFunctorBase.h"
22 #include "Message/Communicate.h"
23 
24 namespace qmcplusplus
25 {
26 /** A numerical functor
27  *
28  * implements interfaces to be used for Jastrow functions and replaces  CubicBsplineSingle.
29  * Template parameters
30  * - RT real data type
31  * - FNOUT final numerical functor
32  *  An example is in CBSOBuilder.h which uses CubicBspline
33  *  typedef CubicBspline<RealType,LINEAR_1DGRID,FIRSTDERIV_CONSTRAINTS> SplineEngineType;
34  *  typedef CubicSplineSingle<RealType,SplineEngineType> RadialOrbitalType;
35  */
36 template<typename RT, typename FNOUT>
37 struct CubicSplineSingle : public OptimizableFunctorBase
38 {
39   ///typedef for the value_type
40   typedef RT value_type;
41   ///typedef of the source functor
42   typedef OptimizableFunctorBase FNIN;
43   ///typedef for the grid
44   typedef OneDimGridBase<real_type> grid_type;
45 
46   // mmorales: until I figure out how to go around this
47   int NumParams;
48   FNIN* InFunc;
49   FNOUT OutFunc;
50   int NumGridPoints;
51   real_type Rmax;
52   real_type GridDelta;
53   real_type Y;
54   real_type dY;
55   real_type d2Y;
56   ///constructor
CubicSplineSingleCubicSplineSingle57   CubicSplineSingle() : InFunc(0) {}
58 
CubicSplineSingleCubicSplineSingle59   CubicSplineSingle(const CubicSplineSingle& old)
60       : OptimizableFunctorBase(old), NumGridPoints(old.NumGridPoints), Rmax(old.Rmax), GridDelta(old.GridDelta)
61   {
62     NumParams = 0;
63     if (old.InFunc)
64     {
65       initialize(old.InFunc->makeClone(), old.Rmax, old.NumGridPoints);
66     }
67   }
68 
makeCloneCubicSplineSingle69   OptimizableFunctorBase* makeClone() const { return new CubicSplineSingle<RT, FNOUT>(*this); }
70 
71   ///constructor with arguments
CubicSplineSingleCubicSplineSingle72   CubicSplineSingle(FNIN* in_, grid_type* agrid) : InFunc(in_) { initialize(in_, agrid); }
73   ///constructor with arguments
CubicSplineSingleCubicSplineSingle74   CubicSplineSingle(FNIN* in_, real_type rc, int npts) : InFunc(in_) { initialize(in_, rc, npts); }
75   ///set the input, analytic function
setInFuncCubicSplineSingle76   void setInFunc(FNIN* in_) { InFunc = in_; }
77 
reportStatusCubicSplineSingle78   void reportStatus(std::ostream& os)
79   {
80     //myVars.print(os);
81   }
82 
isOptimizableCubicSplineSingle83   bool isOptimizable() { return false; }
84 
85   /** evaluate everything: value, first, second and third derivatives
86    */
evaluateCubicSplineSingle87   inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3)
88   {
89     std::cerr << "Third derivative not implemented for CubicSplineSingle.\n";
90     return OutFunc.splint(r, dudr, d2udr2);
91   }
92 
93 
94   /** evaluate everything: value, first and second derivatives
95    */
evaluateCubicSplineSingle96   inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) { return OutFunc.splint(r, dudr, d2udr2); }
97 
98   /** evaluate value only
99    */
evaluateCubicSplineSingle100   inline real_type evaluate(real_type r) { return OutFunc.splint(r); }
101 
102   /** evaluate value only
103    *
104    * Function required for SphericalBasisSet
105    */
evaluateCubicSplineSingle106   inline real_type evaluate(real_type r, real_type rinv) { return Y = OutFunc.splint(r); }
107 
108   /** evaluate everything: value, first and second derivatives
109    *
110    * Function required for SphericalBasisSet
111    */
evaluateAllCubicSplineSingle112   inline real_type evaluateAll(real_type r, real_type rinv) { return Y = OutFunc.splint(r, dY, d2Y); }
113 
114   /** implement the virtual function of OptimizableFunctorBase */
fCubicSplineSingle115   inline real_type f(real_type r) { return OutFunc.splint(r); }
116 
117   /** implement the virtual function of OptimizableFunctorBase  */
dfCubicSplineSingle118   inline real_type df(real_type r)
119   {
120     real_type dudr, d2udr2;
121     OutFunc.splint(r, dudr, d2udr2);
122     return dudr;
123   }
124 
evaluateVCubicSplineSingle125   inline real_type evaluateV(const int iat,
126                              const int iStart,
127                              const int iEnd,
128                              const real_type* restrict _distArray,
129                              real_type* restrict distArrayCompressed) const
130   {
131     // need to actually implement this!
132     return real_type(0);
133   }
134 
evaluateVGLCubicSplineSingle135   inline void evaluateVGL(const int iat,
136                           const int iStart,
137                           const int iEnd,
138                           const real_type* distArray,
139                           real_type* restrict valArray,
140                           real_type* restrict gradArray,
141                           real_type* restrict laplArray,
142                           real_type* restrict distArrayCompressed,
143                           int* restrict distIndices) const
144   {
145     // need to actually implement this!
146   }
147 
putCubicSplineSingle148   bool put(xmlNodePtr cur)
149   {
150     bool s = false;
151     if (InFunc)
152       s = InFunc->put(cur);
153     return s;
154   }
155 
checkInVariablesCubicSplineSingle156   void checkInVariables(opt_variables_type& active)
157   {
158     if (InFunc)
159       InFunc->checkInVariables(active);
160   }
161 
checkOutVariablesCubicSplineSingle162   void checkOutVariables(const opt_variables_type& active)
163   {
164     if (InFunc)
165       InFunc->checkOutVariables(active);
166   }
167 
168   ///reset the input/output function
resetParametersCubicSplineSingle169   void resetParameters(const opt_variables_type& active)
170   {
171     if (InFunc)
172     {
173       InFunc->resetParameters(active);
174       reset();
175     }
176   }
177 
printCubicSplineSingle178   void print(std::ostream& os)
179   {
180     real_type r = 0;
181     for (int i = 0; i < NumGridPoints; i++, r += GridDelta)
182       os << r << " " << OutFunc.splint(r) << std::endl;
183   }
184 
185   ///set the input, analytic function
initializeCubicSplineSingle186   void initialize(FNIN* in_, grid_type* agrid) { initialize(in_, agrid->rmax(), agrid->size()); }
187 
initializeCubicSplineSingle188   void initialize(FNIN* in_, real_type rmax, int npts)
189   {
190     InFunc        = in_;
191     Rmax          = rmax;
192     NumGridPoints = npts;
193     GridDelta     = Rmax / static_cast<real_type>(NumGridPoints - 1);
194     reset();
195   }
196 
resetCubicSplineSingle197   void reset()
198   {
199     if (InFunc)
200     {
201       typename FNOUT::container_type datain(NumGridPoints);
202       real_type r = 0;
203       for (int i = 0; i < NumGridPoints; i++, r += GridDelta)
204       {
205         datain[i] = InFunc->f(r);
206       }
207       OutFunc.Init(0.0, Rmax, datain, true, InFunc->df(0.0), 0.0);
208     }
209     else
210     {
211       APP_ABORT("CubicSplineSingle::reset has no input functor");
212     }
213   }
214 };
215 
216 
217 template<typename RT>
218 struct CubicSplineBasisSet : public OptimizableFunctorBase
219 {
220   ///typedef of the source functor
221   typedef OptimizableFunctorBase FNIN;
222   ///typedef for the argument
223   typedef CubicBspline<RT, LINEAR_1DGRID, FIRSTDERIV_CONSTRAINTS> FNOUT;
224   ///typedef for the grid
225   typedef OneDimGridBase<real_type> grid_type;
226 
227   FNIN* InFunc;
228   FNOUT* OutFunc;
229   int NumGridPoints;
230   real_type Rmax;
231   real_type GridDelta;
232 
233   ///constructor
CubicSplineBasisSetCubicSplineBasisSet234   CubicSplineBasisSet() : InFunc(0), OutFunc(0) {}
235   ///constructor with arguments
CubicSplineBasisSetCubicSplineBasisSet236   CubicSplineBasisSet(FNIN* in_, grid_type* agrid) { initialize(in_, agrid); }
237   ///set the input, analytic function
setInFuncCubicSplineBasisSet238   void setInFunc(FNIN* in_) { InFunc = in_; }
239   ///set the output numerical function
setOutFuncCubicSplineBasisSet240   void setOutFunc(FNOUT* out_) { OutFunc = out_; }
241   ///reset the input/output function
resetParametersCubicSplineBasisSet242   void resetParameters(const opt_variables_type& active)
243   {
244     if (!InFunc)
245       APP_ABORT("CubicSplineBasisSet::resetParameters failed due to null input function ");
246     InFunc->resetParameters(active);
247     reset();
248   }
249 
resetCubicSplineBasisSet250   void reset()
251   {
252     if (!OutFunc)
253       OutFunc = new FNOUT;
254     typename FNOUT::container_type datain(NumGridPoints);
255     real_type r = 0;
256     for (int i = 0; i < NumGridPoints; i++, r += GridDelta)
257       datain[i] = InFunc->f(r);
258     OutFunc->Init(0.0, Rmax, datain, true, InFunc->df(0.0), 0.0);
259   }
260 
261   /** evaluate everything: value, first and second derivatives
262   */
evaluateCubicSplineBasisSet263   inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
264   {
265     return OutFunc->splint(r, dudr, d2udr2);
266   }
267 
268   /** evaluate value only
269   */
evaluateCubicSplineBasisSet270   inline real_type evaluate(real_type r) { return OutFunc->splint(r); }
271 
272   /** implement the virtual function of OptimizableFunctorBase */
fCubicSplineBasisSet273   real_type f(real_type r) { return OutFunc->splint(r); }
274 
275   /** implement the virtual function of OptimizableFunctorBase  */
dfCubicSplineBasisSet276   real_type df(real_type r)
277   {
278     real_type dudr, d2udr2;
279     OutFunc->splint(r, dudr, d2udr2);
280     return dudr;
281   }
282 
putCubicSplineBasisSet283   bool put(xmlNodePtr cur) { return InFunc->put(cur); }
284 
printCubicSplineBasisSet285   void print(std::ostream& os)
286   {
287     real_type r = 0;
288     for (int i = 0; i < NumGridPoints; i++, r += GridDelta)
289       os << r << " " << OutFunc->splint(r) << std::endl;
290   }
291 
292   ///set the input, analytic function
initializeCubicSplineBasisSet293   void initialize(FNIN* in_, grid_type* agrid)
294   {
295     Rmax          = agrid->rmax();
296     NumGridPoints = agrid->size();
297     GridDelta     = Rmax / static_cast<real_type>(NumGridPoints - 1);
298     InFunc        = in_;
299     reset();
300   }
301 };
302 
303 //  /** A numerical functor using cubic bspline
304 //   *
305 //   * This class is replaced by CubicSplineSingle<RT,FNOUT>
306 //   */
307 //  template <typename RT>
308 //    struct CubicBsplineSingle: public OptimizableFunctorBase<RT> {
309 //
310 //      ///typedef of the source functor
311 //      typedef OptimizableFunctorBase<RT> FNIN;
312 //      ///typedef for the argument
313 //      typedef typename FNIN::real_type real_type;
314 //      ///typedef for OptimizableSetType
315 //      typedef typename FNIN::OptimizableSetType OptimizableSetType;
316 //      ///typedef for the argument
317 //      typedef CubicBspline<RT,LINEAR_1DGRID,FIRSTDERIV_CONSTRAINTS> FNOUT;
318 //      ///typedef for the grid
319 //      typedef OneDimGridBase<real_type> grid_type;
320 //
321 //      FNIN *InFunc;
322 //      FNOUT *OutFunc;
323 //      int NumGridPoints;
324 //      real_type Rmax;
325 //      real_type GridDelta;
326 //      real_type Y;
327 //      real_type dY;
328 //      real_type d2Y;
329 //
330 //
331 //      ///constructor
332 //      CubicBsplineSingle(): InFunc(0), OutFunc(0) { }
333 //      ///constructor with arguments
334 //      CubicBsplineSingle(FNIN* in_, grid_type* agrid): InFunc(0), OutFunc(0) {
335 //        initialize(in_,agrid);
336 //      }
337 //      ///constructor with arguments
338 //      CubicBsplineSingle(FNIN* in_, real_type rc, int npts):InFunc(0), OutFunc(0){
339 //        initialize(in_,rc,npts);
340 //      }
341 //      ///set the input, analytic function
342 //      void setInFunc(FNIN* in_) { InFunc=in_;}
343 //      ///set the output numerical function
344 //      void setOutFunc(FNOUT* out_) { OutFunc=out_;}
345 //
346 //      /** evaluate everything: value, first and second derivatives
347 //       */
348 //      inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) {
349 //        return OutFunc->splint(r,dudr,d2udr2);
350 //      }
351 //
352 //      /** evaluate value only
353 //       */
354 //      inline real_type evaluate(real_type r) {
355 //        return OutFunc->splint(r);
356 //      }
357 //
358 //      /** evaluate value only
359 //       *
360 //       * Function required for SphericalBasisSet
361 //       */
362 //      inline real_type evaluate(real_type r, real_type rinv)
363 //      {
364 //        return Y=OutFunc->splint(r);
365 //      }
366 //
367 //      /** evaluate everything: value, first and second derivatives
368 //       *
369 //       * Function required for SphericalBasisSet
370 //       */
371 //      inline real_type evaluateAll(real_type r, real_type rinv)
372 //      {
373 //        return Y=OutFunc->splint(r,dY,d2Y);
374 //      }
375 //
376 //      /** implement the virtual function of OptimizableFunctorBase */
377 //      inline real_type f(real_type r)
378 //      {
379 //        return OutFunc->splint(r);
380 //      }
381 //
382 //      /** implement the virtual function of OptimizableFunctorBase  */
383 //      inline real_type df(real_type r) {
384 //        real_type dudr,d2udr2;
385 //        OutFunc->splint(r,dudr,d2udr2);
386 //        return dudr;
387 //      }
388 //
389 //      bool put(xmlNodePtr cur)
390 //      {
391 //        if(InFunc)
392 //          return InFunc->put(cur);
393 //        else
394 //          return false;
395 //      }
396 //
397 //      void addOptimizables(OptimizableSetType& vlist)
398 //      {
399 //        if(InFunc)
400 //          InFunc->addOptimizables(vlist);
401 //      }
402 //
403 //      ///reset the input/output function
404 //      void resetParameters(OptimizableSetType& optVariables)
405 //      {
406 //        if(!InFunc)
407 //        {
408 //          app_error() << "  CubicSplineJastrow::reset failed due to null input function " << std::endl;
409 //          OHMMS::Controller->abort();
410 //        }
411 //        InFunc->resetParameters(optVariables);
412 //        resetInternals();
413 //      }
414 //
415 //      void print(std::ostream& os) {
416 //        real_type r=0;
417 //        for(int i=0; i<NumGridPoints; i++, r+=GridDelta)
418 //          os << r << " " << OutFunc->splint(r) << std::endl;
419 //      }
420 //
421 //      ///set the input, analytic function
422 //      void initialize(FNIN* in_, grid_type* agrid) {
423 //        initialize(in_,agrid->rmax(),agrid->size());
424 //      }
425 //      void initialize(FNIN* in_, real_type rmax, int npts)
426 //      {
427 //        InFunc=in_;
428 //        Rmax=rmax;
429 //        NumGridPoints=npts;
430 //        GridDelta=Rmax/static_cast<real_type>(NumGridPoints-1);
431 //        resetInternals();
432 //      }
433 //
434 //      void resetInternals()
435 //      {
436 //        if(!OutFunc) OutFunc = new FNOUT;
437 //        typename FNOUT::container_type datain(NumGridPoints);
438 //        real_type r=0;
439 //        for(int i=0; i<NumGridPoints; i++, r+=GridDelta)
440 //        {
441 //          datain[i] = InFunc->f(r);
442 //        }
443 //        OutFunc->Init(0.0,Rmax,datain,true,InFunc->df(0.0),0.0);
444 //      }
445 //
446 //    };
447 //
448 //  /** A numerical functor
449 //   *
450 //   * implements interfaces to be used for Jastrow functions
451 //   * - OneBodyJastrow<NumericalJastrow>
452 //   * - TwoBodyJastrow<NumericalJastrow>
453 //   */
454 //  template <class RT>
455 //    struct SplineJastrow: public OptimizableFunctorBase<RT> {
456 //
457 //      ///typedef of the source functor
458 //      typedef OptimizableFunctorBase<RT> FNIN;
459 //      ///typedef for the argument
460 //      typedef typename FNIN::real_type real_type;
461 //      ///typedef of the target functor
462 //      typedef OneDimCubicSpline<real_type,real_type>  FNOUT;
463 //
464 //      real_type Rmax;
465 //      FNIN *InFunc;
466 //      FNOUT *OutFunc;
467 //
468 //      ///constructor
469 //      SplineJastrow(): InFunc(0), OutFunc(0) { }
470 //      ///constructor with arguments
471 //      SplineJastrow(FNIN* in_, typename FNOUT::grid_type* agrid){
472 //        initialize(in_,agrid);
473 //      }
474 //      ///set the input, analytic function
475 //      void setInFunc(FNIN* in_) { InFunc=in_;}
476 //      ///set the output numerical function
477 //      void setOutFunc(FNOUT* out_) { OutFunc=out_;}
478 //      ///reset the input/output function
479 //      inline void reset() {
480 //        InFunc->reset();
481 //        //reference to the output functions grid
482 //        const typename FNOUT::grid_type& grid = OutFunc->grid();
483 //        //set cutoff function
484 //        int last=grid.size()-1;
485 //        for(int i=0; i<grid.size(); i++) {
486 //          (*OutFunc)(i) = InFunc->f(grid(i));
487 ////	  std::cout << grid(i) << "   " << (*OutFunc)(i) << std::endl;
488 //        }
489 //	(*OutFunc)(last)=0.0;
490 //        //boundary conditions
491 //        real_type deriv1=InFunc->df(grid(0));
492 //        real_type deriv2=0.0;
493 //        OutFunc->spline(0,deriv1,last,deriv2);
494 //        Rmax=grid(last);
495 //      }
496 //
497 //      /** evaluate everything: value, first and second derivatives
498 //      */
499 //      inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) {
500 //        return OutFunc->splint(r,dudr,d2udr2);
501 //      }
502 //
503 //      /** evaluate value only
504 //      */
505 //      inline real_type evaluate(real_type r) {
506 //        return OutFunc->splint(r);
507 //      }
508 //
509 //      /** implement the virtual function of OptimizableFunctorBase */
510 //      real_type f(real_type r) {
511 //        return OutFunc->splint(r);
512 //      }
513 //
514 //      /** implement the virtual function of OptimizableFunctorBase  */
515 //      real_type df(real_type r) {
516 //        real_type dudr,d2udr2;
517 //        OutFunc->splint(r,dudr,d2udr2);
518 //        return dudr;
519 //      }
520 //
521 //      bool put(xmlNodePtr cur)
522 //      {
523 //        return InFunc->put(cur);
524 //      }
525 //
526 //      void addOptimizables( VarRegistry<real_type>& vlist)
527 //      {
528 //        InFunc->addOptimizables(vlist);
529 //      }
530 //      void print(std::ostream& os) {
531 //        const typename FNOUT::grid_type& grid = OutFunc->grid();
532 //        for(int i=0; i<grid.size(); i++) {
533 //          std::cout << grid(i) << " " << (*OutFunc)(i) << std::endl;
534 //        }
535 //      }
536 //
537 //      ///set the input, analytic function
538 //      void initialize(FNIN* in_, typename FNOUT::grid_type* agrid) {
539 //        InFunc=in_;
540 //        setOutFunc(new FNOUT(agrid));
541 //        reset();
542 //      }
543 //    };
544 //
545 } // namespace qmcplusplus
546 #endif
547