Inheritance diagram for PolynomialFunctor3D:

Collaboration diagram for PolynomialFunctor3D:

Public Types
using	value_type = real_type

Public Types inherited from OptimizableFunctorBase
using	real_type = optimize::VariableSet::real_type
	typedef for real values More...

using	opt_variables_type = optimize::VariableSet
	typedef for variableset: this is going to be replaced More...

Public Member Functions
	PolynomialFunctor3D (const std::string &my_name, real_type ee_cusp=0.0, real_type eI_cusp=0.0)
	constructor More...

OptimizableFunctorBase *	makeClone () const override
	create a clone of this object More...

void	resize (int neI, int nee)

void	reset () override
	reset function More...

void	reset_gamma ()

real_type	evaluate (real_type r_12, real_type r_1I, real_type r_2I) const

real_type	evaluateV (int Nptcl, const real_type restrict r_12_array, const real_type restrict r_1I_array, const real_type *restrict r_2I_array) const

real_type	evaluate (real_type r_12, real_type r_1I, real_type r_2I, TinyVector< real_type, 3 > &grad, Tensor< real_type, 3 > &hess) const

void	evaluateVGL (int Nptcl, const real_type restrict r_12_array, const real_type restrict r_1I_array, const real_type restrict r_2I_array, real_type restrict val_array, real_type restrict grad0_array, real_type restrict grad1_array, real_type restrict grad2_array, real_type restrict hess00_array, real_type restrict hess11_array, real_type restrict hess22_array, real_type restrict hess01_array, real_type restrict hess02_array) const

real_type	evaluate (const real_type r_12, const real_type r_1I, const real_type r_2I, TinyVector< real_type, 3 > &grad, Tensor< real_type, 3 > &hess, TinyVector< Tensor< real_type, 3 >, 3 > &d3)

real_type	evaluate (const real_type r, const real_type rinv)

bool	evaluateDerivativesFD (const real_type r_12, const real_type r_1I, const real_type r_2I, std::vector< double > &d_vals, std::vector< TinyVector< real_type, 3 >> &d_grads, std::vector< Tensor< real_type, 3 >> &d_hess)

bool	evaluateDerivatives (const real_type r_12, const real_type r_1I, const real_type r_2I, std::vector< real_type > &d_vals)
	calculate derivatives with respect to polynomial parameters More...

bool	evaluateDerivatives (const real_type r_12, const real_type r_1I, const real_type r_2I, std::vector< real_type > &d_vals, std::vector< TinyVector< real_type, 3 >> &d_grads, std::vector< Tensor< real_type, 3 >> &d_hess)

real_type	f (real_type r) override

real_type	df (real_type r) override

bool	put (xmlNodePtr cur) override
	process xmlnode and registers variables to optimize More...

void	resetParametersExclusive (const opt_variables_type &active) override
	reset the parameters during optimizations. More...

void	checkInVariablesExclusive (opt_variables_type &active) override
	check in variational parameters to the global list of parameters used by the optimizer. More...

void	checkOutVariables (const opt_variables_type &active) override
	check out variational optimizable variables More...

void	print (std::ostream &os)

int	getNumParameters ()

Public Member Functions inherited from OptimizableFunctorBase
	OptimizableFunctorBase (const std::string &name="")
	default constructor More...

virtual	~OptimizableFunctorBase ()=default
	virtual destrutor More...

void	getIndex (const opt_variables_type &active)

virtual real_type	f (real_type r)=0
	evaluate the value at r More...

virtual real_type	df (real_type r)=0
	evaluate the first derivative More...

virtual void	setDensity (real_type n)
	empty virtual function to help builder classes More...

virtual void	setCusp (real_type cusp)
	empty virtual function to help builder classes More...

virtual void	setPeriodic (bool periodic)
	empty virtual function to help builder classes More...

virtual bool	evaluateDerivatives (real_type r, std::vector< qmcplusplus::TinyVector< real_type, 3 >> &derivs)

virtual bool	evaluateDerivatives (real_type r, std::vector< real_type > &derivs)

virtual void	setGridManager (bool willmanage)

virtual void	checkInVariablesExclusive (opt_variables_type &active)=0
	check in variational parameters to the global list of parameters used by the optimizer. More...

virtual void	resetParametersExclusive (const opt_variables_type &active)=0
	reset the parameters during optimizations More...

Public Member Functions inherited from OptimizableObject
	OptimizableObject (const std::string &name)

const std::string &	getName () const

bool	isOptimized () const

virtual void	reportStatus (std::ostream &os)
	print the state, e.g., optimizables More...

void	setOptimization (bool state)

virtual void	writeVariationalParameters (hdf_archive &hout)
	Write the variational parameters for this object to the VP HDF file. More...

virtual void	readVariationalParameters (hdf_archive &hin)
	Read the variational parameters for this object from the VP HDF file. More...

Public Attributes
int	N_eI

int	N_ee

Array< real_type, 3 >	gamma

std::vector< int >	GammaPerm

Array< int, 3 >	index

std::vector< bool >	IndepVar

std::vector< real_type >	GammaVec

std::vector< real_type >	dval_Vec

std::vector< TinyVector< real_type, 3 > >	dgrad_Vec

std::vector< Tensor< real_type, 3 > >	dhess_Vec

int	NumConstraints

int	NumGamma

Matrix< real_type >	ConstraintMatrix

std::vector< real_type >	Parameters

std::vector< real_type >	d_valsFD

std::vector< TinyVector< real_type, 3 > >	d_gradsFD

std::vector< Tensor< real_type, 3 > >	d_hessFD

std::vector< std::string >	ParameterNames

std::string	iSpecies

std::string	eSpecies1

std::string	eSpecies2

int	ResetCount

const int	C

real_type	scale

bool	notOpt

Public Attributes inherited from OptimizableFunctorBase
real_type	cutoff_radius = 0.0
	maximum cutoff More...

opt_variables_type	myVars
	set of variables to be optimized More...

Detailed Description

Definition at line 30 of file PolynomialFunctor3D.h.

Member Typedef Documentation

◆ value_type

using value_type = real_type

Definition at line 32 of file PolynomialFunctor3D.h.

Constructor & Destructor Documentation

◆ PolynomialFunctor3D()

PolynomialFunctor3D	(	const std::string &	my_name,
		real_type	ee_cusp = `0.0`,
		real_type	eI_cusp = `0.0`
	)

inline

constructor

Definition at line 58 of file PolynomialFunctor3D.h.

References qmcplusplus::abs(), qmcplusplus::app_error(), and OptimizableFunctorBase::cutoff_radius.

Referenced by PolynomialFunctor3D::makeClone().

       : OptimizableFunctorBase(my_name), N_eI(0), N_ee(0), ResetCount(0), C(3), scale(1.0), notOpt(false)
   {
     if (std::abs(ee_cusp) > 0.0 || std::abs(eI_cusp) > 0.0)
     {
       app_error() << "PolynomialFunctor3D does not support nonzero cusp.\n";
       abort();
     }
     cutoff_radius = 0.0;
   }

Member Function Documentation

◆ checkInVariablesExclusive()

void checkInVariablesExclusive ( opt_variables_type & active )

inlineoverridevirtual

check in variational parameters to the global list of parameters used by the optimizer.

Parameters

active a super set of optimizable variables

The existing checkInVariables implementation in WFC/SPO/.. are inclusive and it calls checkInVariables of its members class A: public SPOSet {} class B: public WFC { A objA; checkInVariables() { objA.checkInVariables(); } };

With OptimizableObject, class A: public OptimizableObject {} class B: public OptimizableObject { A objA; checkInVariablesExclusive() { // should not call objA.checkInVariablesExclusive() if objA has been extracted; } }; A vector of OptimizableObject, will be created by calling extractOptimizableObjects(). All the checkInVariablesExclusive() will be called through this vector and thus checkInVariablesExclusive implementation should only handle non-OptimizableObject members.

Implements OptimizableObject.

Definition at line 1073 of file PolynomialFunctor3D.h.

References VariableSet::insertFrom(), OptimizableFunctorBase::myVars, PolynomialFunctor3D::notOpt, and VariableSet::setIndexDefault().

   {
     if (notOpt)
       return;
 
     myVars.setIndexDefault();
     active.insertFrom(myVars);
   }

◆ checkOutVariables()

void checkOutVariables ( const opt_variables_type & active )

inlineoverridevirtual

check out variational optimizable variables

Parameters

active a super set of optimizable variables

Implements OptimizableFunctorBase.

Definition at line 1082 of file PolynomialFunctor3D.h.

References VariableSet::getIndex(), OptimizableFunctorBase::myVars, and PolynomialFunctor3D::notOpt.

   {
     if (notOpt)
       return;
     myVars.getIndex(active);
   }

◆ df()

real_type df ( real_type r )

inlineoverride

Definition at line 969 of file PolynomialFunctor3D.h.

969 { return 0.0; }

◆ evaluate() [1/4]

real_type evaluate	(	real_type	r_12,
		real_type	r_1I,
		real_type	r_2I
	)		const

inline

Definition at line 320 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::C, BLAS::cone, OptimizableFunctorBase::cutoff_radius, BLAS::czero, PolynomialFunctor3D::gamma, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, and PolynomialFunctor3D::N_eI.

Referenced by PolynomialFunctor3D::evaluateDerivativesFD(), and qmcplusplus::TEST_CASE().

   {
     constexpr real_type czero(0);
     constexpr real_type cone(1);
     constexpr real_type chalf(0.5);
 
     const real_type L = chalf * cutoff_radius;
     if (r_1I >= L || r_2I >= L)
       return czero;
     real_type val = czero;
     real_type r2l(cone);
     for (int l = 0; l <= N_eI; l++)
     {
       real_type r2m(r2l);
       for (int m = 0; m <= N_eI; m++)
       {
         real_type r2n(r2m);
         for (int n = 0; n <= N_ee; n++)
         {
           val += gamma(l, m, n) * r2n;
           r2n *= r_12;
         }
         r2m *= r_2I;
       }
       r2l *= r_1I;
     }
     for (int i = 0; i < C; i++)
       val *= (r_1I - L) * (r_2I - L);
     return val;
   }

◆ evaluate() [2/4]

real_type evaluate	(	real_type	r_12,
		real_type	r_1I,
		real_type	r_2I,
		TinyVector< real_type, 3 > &	grad,
		Tensor< real_type, 3 > &	hess
	)		const

inline

Definition at line 396 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::C, BLAS::cone, OptimizableFunctorBase::cutoff_radius, BLAS::czero, PolynomialFunctor3D::gamma, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, and PolynomialFunctor3D::N_eI.

   {
     constexpr real_type czero(0);
     constexpr real_type cone(1);
     constexpr real_type chalf(0.5);
     constexpr real_type ctwo(2);
 
     const real_type L = chalf * cutoff_radius;
     if (r_1I >= L || r_2I >= L)
     {
       grad = czero;
       hess = czero;
       return czero;
     }
     real_type val = czero;
     grad          = czero;
     hess          = czero;
     real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero);
     for (int l = 0; l <= N_eI; l++)
     {
       real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero);
       for (int m = 0; m <= N_eI; m++)
       {
         real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero);
         for (int n = 0; n <= N_ee; n++)
         {
           const real_type g    = gamma(l, m, n);
           const real_type g00x = g * r2l * r2m;
           const real_type g10x = g * r2l_1 * r2m;
           const real_type g01x = g * r2l * r2m_1;
           const real_type gxx0 = g * r2n;
 
           val += g00x * r2n;
           grad[0] += g00x * r2n_1;
           grad[1] += g10x * r2n;
           grad[2] += g01x * r2n;
           hess(0, 0) += g00x * r2n_2;
           hess(0, 1) += g10x * r2n_1;
           hess(0, 2) += g01x * r2n_1;
           hess(1, 1) += gxx0 * r2l_2 * r2m;
           hess(1, 2) += gxx0 * r2l_1 * r2m_1;
           hess(2, 2) += gxx0 * r2l * r2m_2;
           nf += cone;
           r2n_2 = r2n_1 * nf;
           r2n_1 = r2n * nf;
           r2n *= r_12;
         }
         mf += cone;
         r2m_2 = r2m_1 * mf;
         r2m_1 = r2m * mf;
         r2m *= r_2I;
       }
       lf += cone;
       r2l_2 = r2l_1 * lf;
       r2l_1 = r2l * lf;
       r2l *= r_1I;
     }
 
     const real_type r_2I_minus_L = r_2I - L;
     const real_type r_1I_minus_L = r_1I - L;
     const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L;
     for (int i = 0; i < C; i++)
     {
       hess(0, 0) = both_minus_L * hess(0, 0);
       hess(0, 1) = both_minus_L * hess(0, 1) + r_2I_minus_L * grad[0];
       hess(0, 2) = both_minus_L * hess(0, 2) + r_1I_minus_L * grad[0];
       hess(1, 1) = both_minus_L * hess(1, 1) + ctwo * r_2I_minus_L * grad[1];
       hess(1, 2) = both_minus_L * hess(1, 2) + r_1I_minus_L * grad[1] + r_2I_minus_L * grad[2] + val;
       hess(2, 2) = both_minus_L * hess(2, 2) + ctwo * r_1I_minus_L * grad[2];
       grad[0]    = both_minus_L * grad[0];
       grad[1]    = both_minus_L * grad[1] + r_2I_minus_L * val;
       grad[2]    = both_minus_L * grad[2] + r_1I_minus_L * val;
       val *= both_minus_L;
     }
     hess(1, 0) = hess(0, 1);
     hess(2, 0) = hess(0, 2);
     hess(2, 1) = hess(1, 2);
     return val;
   }

◆ evaluate() [3/4]

real_type evaluate	(	const real_type	r_12,
		const real_type	r_1I,
		const real_type	r_2I,
		TinyVector< real_type, 3 > &	grad,
		Tensor< real_type, 3 > &	hess,
		TinyVector< Tensor< real_type, 3 >, 3 > &	d3
	)

inline

Definition at line 589 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::C, OptimizableFunctorBase::cutoff_radius, PolynomialFunctor3D::gamma, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, and PolynomialFunctor3D::N_eI.

   {
     grad              = 0.0;
     hess              = 0.0;
     d3                = grad;
     const real_type L = 0.5 * cutoff_radius;
     if (r_1I >= L || r_2I >= L)
       return 0.0;
     real_type val = 0.0;
     real_type r2l(1.0), r2l_1(0.0), r2l_2(0.0), r2l_3(0.0), lf(0.0);
     for (int l = 0; l <= N_eI; l++)
     {
       real_type r2m(1.0), r2m_1(0.0), r2m_2(0.0), r2m_3(0.0), mf(0.0);
       for (int m = 0; m <= N_eI; m++)
       {
         real_type r2n(1.0), r2n_1(0.0), r2n_2(0.0), r2n_3(0.0), nf(0.0);
         for (int n = 0; n <= N_ee; n++)
         {
           real_type g = gamma(l, m, n);
           val += g * r2l * r2m * r2n;
           grad[0] += nf * g * r2l * r2m * r2n_1;
           grad[1] += lf * g * r2l_1 * r2m * r2n;
           grad[2] += mf * g * r2l * r2m_1 * r2n;
           hess(0, 0) += nf * (nf - 1.0) * g * r2l * r2m * r2n_2;
           hess(0, 1) += nf * lf * g * r2l_1 * r2m * r2n_1;
           hess(0, 2) += nf * mf * g * r2l * r2m_1 * r2n_1;
           hess(1, 1) += lf * (lf - 1.0) * g * r2l_2 * r2m * r2n;
           hess(1, 2) += lf * mf * g * r2l_1 * r2m_1 * r2n;
           hess(2, 2) += mf * (mf - 1.0) * g * r2l * r2m_2 * r2n;
           d3[0](0, 0) += nf * (nf - 1.0) * (nf - 2.0) * g * r2l * r2m * r2n_3;
           d3[0](0, 1) += nf * (nf - 1.0) * lf * g * r2l_1 * r2m * r2n_2;
           d3[0](0, 2) += nf * (nf - 1.0) * mf * g * r2l * r2m_1 * r2n_2;
           d3[0](1, 1) += nf * lf * (lf - 1.0) * g * r2l_2 * r2m * r2n_1;
           d3[0](1, 2) += nf * lf * mf * g * r2l_1 * r2m_1 * r2n_1;
           d3[0](2, 2) += nf * mf * (mf - 1.0) * g * r2l * r2m_2 * r2n_1;
           d3[1](1, 1) += lf * (lf - 1.0) * (lf - 2.0) * g * r2l_3 * r2m * r2n;
           d3[1](1, 2) += lf * (lf - 1.0) * mf * g * r2l_2 * r2m_1 * r2n;
           d3[1](2, 2) += lf * mf * (mf - 1.0) * g * r2l_1 * r2m_2 * r2n;
           d3[2](2, 2) += mf * (mf - 1.0) * (mf - 2.0) * g * r2l * r2m_3 * r2n;
           r2n_3 = r2n_2;
           r2n_2 = r2n_1;
           r2n_1 = r2n;
           r2n *= r_12;
           nf += 1.0;
         }
         r2m_3 = r2m_2;
         r2m_2 = r2m_1;
         r2m_1 = r2m;
         r2m *= r_2I;
         mf += 1.0;
       }
       r2l_3 = r2l_2;
       r2l_2 = r2l_1;
       r2l_1 = r2l;
       r2l *= r_1I;
       lf += 1.0;
     }
     for (int i = 0; i < C; i++)
     {
       d3[0](0, 0) = (r_1I - L) * (r_2I - L) * d3[0](0, 0);
       d3[0](0, 1) = (r_1I - L) * (r_2I - L) * d3[0](0, 1) + (r_2I - L) * hess(0, 0);
       d3[0](0, 2) = (r_1I - L) * (r_2I - L) * d3[0](0, 2) + (r_1I - L) * hess(0, 0);
       d3[0](1, 1) = (r_1I - L) * (r_2I - L) * d3[0](1, 1) + 2.0 * (r_2I - L) * hess(0, 1);
       d3[0](1, 2) = (r_1I - L) * (r_2I - L) * d3[0](1, 2) + (r_1I - L) * hess(0, 1) + (r_2I - L) * hess(0, 2) + grad[0];
       d3[0](2, 2) = (r_1I - L) * (r_2I - L) * d3[0](2, 2) + 2.0 * (r_1I - L) * hess(0, 2);
       d3[1](1, 1) = (r_1I - L) * (r_2I - L) * d3[1](1, 1) + 3.0 * (r_2I - L) * hess(1, 1);
       d3[1](1, 2) = (r_1I - L) * (r_2I - L) * d3[1](1, 2) + 2.0 * (r_2I - L) * hess(1, 2) + 2.0 * grad[1] +
           (r_1I - L) * hess(1, 1);
       d3[1](2, 2) = (r_1I - L) * (r_2I - L) * d3[1](2, 2) + 2.0 * (r_1I - L) * hess(1, 2) + 2.0 * grad[2] +
           (r_2I - L) * hess(2, 2);
       d3[2](2, 2) = (r_1I - L) * (r_2I - L) * d3[2](2, 2) + 3.0 * (r_1I - L) * hess(2, 2);
       hess(0, 0)  = (r_1I - L) * (r_2I - L) * hess(0, 0);
       hess(0, 1)  = (r_1I - L) * (r_2I - L) * hess(0, 1) + (r_2I - L) * grad[0];
       hess(0, 2)  = (r_1I - L) * (r_2I - L) * hess(0, 2) + (r_1I - L) * grad[0];
       hess(1, 1)  = (r_1I - L) * (r_2I - L) * hess(1, 1) + 2.0 * (r_2I - L) * grad[1];
       hess(1, 2)  = (r_1I - L) * (r_2I - L) * hess(1, 2) + (r_1I - L) * grad[1] + (r_2I - L) * grad[2] + val;
       hess(2, 2)  = (r_1I - L) * (r_2I - L) * hess(2, 2) + 2.0 * (r_1I - L) * grad[2];
       grad[0]     = (r_1I - L) * (r_2I - L) * grad[0];
       grad[1]     = (r_1I - L) * (r_2I - L) * grad[1] + (r_2I - L) * val;
       grad[2]     = (r_1I - L) * (r_2I - L) * grad[2] + (r_1I - L) * val;
       val *= (r_1I - L) * (r_2I - L);
     }
     hess(1, 0)  = hess(0, 1);
     hess(2, 0)  = hess(0, 2);
     hess(2, 1)  = hess(1, 2);
     d3[0](1, 0) = d3[0](0, 1);
     d3[0](2, 0) = d3[0](0, 2);
     d3[0](2, 1) = d3[0](1, 2);
     d3[1](0, 0) = d3[0](1, 1);
     d3[1](0, 1) = d3[0](0, 1);
     d3[1](0, 2) = d3[0](1, 2);
     d3[1](1, 0) = d3[0](0, 1);
     d3[1](2, 0) = d3[0](1, 2);
     d3[1](2, 1) = d3[1](1, 2);
     d3[2](0, 0) = d3[0](0, 2);
     d3[2](0, 1) = d3[0](1, 2);
     d3[2](0, 2) = d3[0](2, 2);
     d3[2](1, 0) = d3[0](1, 2);
     d3[2](1, 1) = d3[1](1, 2);
     d3[2](1, 2) = d3[1](2, 2);
     d3[2](2, 0) = d3[0](2, 2);
     d3[2](2, 1) = d3[1](2, 2);
     return val;
   }

◆ evaluate() [4/4]

real_type evaluate	(	const real_type	r,
		const real_type	rinv
	)

inline

Definition at line 700 of file PolynomialFunctor3D.h.

700 { return 0.0; }

◆ evaluateDerivatives() [1/2]

bool evaluateDerivatives	(	const real_type	r_12,
		const real_type	r_1I,
		const real_type	r_2I,
		std::vector< real_type > &	d_vals
	)

inline

calculate derivatives with respect to polynomial parameters

Definition at line 737 of file PolynomialFunctor3D.h.

References qmcplusplus::abs(), PolynomialFunctor3D::C, BLAS::cone, PolynomialFunctor3D::ConstraintMatrix, OptimizableFunctorBase::cutoff_radius, BLAS::czero, PolynomialFunctor3D::dval_Vec, qmcplusplus::Units::charge::e, PolynomialFunctor3D::IndepVar, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, PolynomialFunctor3D::N_eI, PolynomialFunctor3D::NumGamma, and PolynomialFunctor3D::scale.

   {
     const real_type L = 0.5 * cutoff_radius;
     if (r_1I >= L || r_2I >= L)
       return false;
 
     constexpr real_type czero(0);
     constexpr real_type cone(1);
 
     real_type dval_dgamma;
 
     for (int i = 0; i < dval_Vec.size(); i++)
       dval_Vec[i]  = czero;
 
     const real_type r_2I_minus_L = r_2I - L;
     const real_type r_1I_minus_L = r_1I - L;
     const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L;
 
     real_type r2l(cone);
     for (int l = 0; l <= N_eI; l++)
     {
       real_type r2m(cone);
       for (int m = 0; m <= N_eI; m++)
       {
         int num;
         if (m > l)
           num = ((2 * N_eI - l + 3) * l / 2 + m - l) * (N_ee + 1);
         else
           num = ((2 * N_eI - m + 3) * m / 2 + l - m) * (N_ee + 1);
         real_type r2n(cone);
         for (int n = 0; n <= N_ee; n++, num++)
         {
           dval_dgamma        = r2l * r2m * r2n;
           for (int i = 0; i < C; i++)
             dval_dgamma *= both_minus_L;
 
           // Now, pack into vectors
           dval_Vec[num] += scale * dval_dgamma;
           r2n *= r_12;
         }
         r2m *= r_2I;
       }
       r2l *= r_1I;
     }
     // for (int i=0; i<dval_Vec.size(); i++)
     //  fprintf (stderr, "dval_Vec[%d] = %12.6e\n", i, dval_Vec[i]);
     ///////////////////////////////////////////
     // Now, compensate for constraint matrix //
     ///////////////////////////////////////////
     std::fill(d_vals.begin(), d_vals.end(), 0.0);
     int var = 0;
     for (int i = 0; i < NumGamma; i++)
       if (IndepVar[i])
       {
         d_vals[var]  = dval_Vec[i];
         var++;
       }
     int constraint = 0;
     for (int i = 0; i < NumGamma; i++)
     {
       if (!IndepVar[i])
       {
         int indep_var = 0;
         for (int j = 0; j < NumGamma; j++)
           if (IndepVar[j])
           {
             d_vals[indep_var] -= ConstraintMatrix(constraint, j) * dval_Vec[i];
             indep_var++;
           }
           else if (i != j)
             assert(std::abs(ConstraintMatrix(constraint, j)) < 1.0e-10);
         constraint++;
       }
     }
     return true;
   }

◆ evaluateDerivatives() [2/2]

bool evaluateDerivatives	(	const real_type	r_12,
		const real_type	r_1I,
		const real_type	r_2I,
		std::vector< real_type > &	d_vals,
		std::vector< TinyVector< real_type, 3 >> &	d_grads,
		std::vector< Tensor< real_type, 3 >> &	d_hess
	)

inline

Definition at line 817 of file PolynomialFunctor3D.h.

   {
     const real_type L = 0.5 * cutoff_radius;
     if (r_1I >= L || r_2I >= L)
       return false;
 
     constexpr real_type czero(0);
     constexpr real_type cone(1);
     constexpr real_type ctwo(2);
 
     real_type dval_dgamma;
     TinyVector<real_type, 3> dgrad_dgamma;
     Tensor<real_type, 3> dhess_dgamma;
 
     for (int i = 0; i < dval_Vec.size(); i++)
     {
       dval_Vec[i]  = czero;
       dgrad_Vec[i] = czero;
       dhess_Vec[i] = czero;
     }
 
     const real_type r_2I_minus_L = r_2I - L;
     const real_type r_1I_minus_L = r_1I - L;
     const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L;
 
     real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero);
     for (int l = 0; l <= N_eI; l++)
     {
       real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero);
       for (int m = 0; m <= N_eI; m++)
       {
         int num;
         if (m > l)
           num = ((2 * N_eI - l + 3) * l / 2 + m - l) * (N_ee + 1);
         else
           num = ((2 * N_eI - m + 3) * m / 2 + l - m) * (N_ee + 1);
         real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero);
         for (int n = 0; n <= N_ee; n++, num++)
         {
           dval_dgamma        = r2l * r2m * r2n;
           dgrad_dgamma[0]    = r2l * r2m * r2n_1;
           dgrad_dgamma[1]    = r2l_1 * r2m * r2n;
           dgrad_dgamma[2]    = r2l * r2m_1 * r2n;
           dhess_dgamma(0, 0) = r2l * r2m * r2n_2;
           dhess_dgamma(0, 1) = r2l_1 * r2m * r2n_1;
           dhess_dgamma(0, 2) = r2l * r2m_1 * r2n_1;
           dhess_dgamma(1, 1) = r2l_2 * r2m * r2n;
           dhess_dgamma(1, 2) = r2l_1 * r2m_1 * r2n;
           dhess_dgamma(2, 2) = r2l * r2m_2 * r2n;
 
           for (int i = 0; i < C; i++)
           {
             dhess_dgamma(0, 0) = both_minus_L * dhess_dgamma(0, 0);
             dhess_dgamma(0, 1) = both_minus_L * dhess_dgamma(0, 1) + r_2I_minus_L * dgrad_dgamma[0];
             dhess_dgamma(0, 2) = both_minus_L * dhess_dgamma(0, 2) + r_1I_minus_L * dgrad_dgamma[0];
             dhess_dgamma(1, 1) = both_minus_L * dhess_dgamma(1, 1) + ctwo * r_2I_minus_L * dgrad_dgamma[1];
             dhess_dgamma(1, 2) = both_minus_L * dhess_dgamma(1, 2) + r_1I_minus_L * dgrad_dgamma[1] +
                 r_2I_minus_L * dgrad_dgamma[2] + dval_dgamma;
             dhess_dgamma(2, 2) = both_minus_L * dhess_dgamma(2, 2) + ctwo * r_1I_minus_L * dgrad_dgamma[2];
             dgrad_dgamma[0]    = both_minus_L * dgrad_dgamma[0];
             dgrad_dgamma[1]    = both_minus_L * dgrad_dgamma[1] + r_2I_minus_L * dval_dgamma;
             dgrad_dgamma[2]    = both_minus_L * dgrad_dgamma[2] + r_1I_minus_L * dval_dgamma;
             dval_dgamma *= both_minus_L;
           }
 
           // Now, pack into vectors
           dval_Vec[num] += scale * dval_dgamma;
           for (int i = 0; i < 3; i++)
           {
             dgrad_Vec[num][i] += scale * dgrad_dgamma[i];
             dhess_Vec[num](i, i) += scale * dhess_dgamma(i, i);
             for (int j = i + 1; j < 3; j++)
             {
               dhess_Vec[num](i, j) += scale * dhess_dgamma(i, j);
               dhess_Vec[num](j, i) = dhess_Vec[num](i, j);
             }
           }
 
           nf += cone;
           r2n_2 = r2n_1 * nf;
           r2n_1 = r2n * nf;
           r2n *= r_12;
         }
         mf += cone;
         r2m_2 = r2m_1 * mf;
         r2m_1 = r2m * mf;
         r2m *= r_2I;
       }
       lf += cone;
       r2l_2 = r2l_1 * lf;
       r2l_1 = r2l * lf;
       r2l *= r_1I;
     }
     // for (int i=0; i<dval_Vec.size(); i++)
     //  fprintf (stderr, "dval_Vec[%d] = %12.6e\n", i, dval_Vec[i]);
     ///////////////////////////////////////////
     // Now, compensate for constraint matrix //
     ///////////////////////////////////////////
     std::fill(d_vals.begin(), d_vals.end(), 0.0);
     int var = 0;
     for (int i = 0; i < NumGamma; i++)
       if (IndepVar[i])
       {
         d_vals[var]  = dval_Vec[i];
         d_grads[var] = dgrad_Vec[i];
         d_hess[var]  = dhess_Vec[i];
         var++;
       }
     int constraint = 0;
     for (int i = 0; i < NumGamma; i++)
     {
       if (!IndepVar[i])
       {
         int indep_var = 0;
         for (int j = 0; j < NumGamma; j++)
           if (IndepVar[j])
           {
             d_vals[indep_var] -= ConstraintMatrix(constraint, j) * dval_Vec[i];
             d_grads[indep_var] -= ConstraintMatrix(constraint, j) * dgrad_Vec[i];
             d_hess[indep_var] -= ConstraintMatrix(constraint, j) * dhess_Vec[i];
             indep_var++;
           }
           else if (i != j)
             assert(std::abs(ConstraintMatrix(constraint, j)) < 1.0e-10);
         constraint++;
       }
     }
     return true;
 #ifdef DEBUG_DERIVS
     evaluateDerivativesFD(r_12, r_1I, r_2I, d_valsFD, d_gradsFD, d_hessFD);
     fprintf(stderr, "Param   Analytic   Finite diffference\n");
     for (int ip = 0; ip < Parameters.size(); ip++)
       fprintf(stderr, "  %3d  %12.6e  %12.6e\n", ip, d_vals[ip], d_valsFD[ip]);
     fprintf(stderr, "Param   Analytic   Finite diffference\n");
     for (int ip = 0; ip < Parameters.size(); ip++)
       fprintf(stderr, "  %3d  %12.6e %12.6e   %12.6e %12.6e   %12.6e %12.6e\n", ip, d_grads[ip][0], d_gradsFD[ip][0],
               d_grads[ip][1], d_gradsFD[ip][1], d_grads[ip][2], d_gradsFD[ip][2]);
     fprintf(stderr, "Param   Analytic   Finite diffference\n");
     for (int ip = 0; ip < Parameters.size(); ip++)
       for (int dim = 0; dim < 3; dim++)
         fprintf(stderr, "  %3d  %12.6e %12.6e   %12.6e %12.6e   %12.6e %12.6e\n", ip, d_hess[ip](0, dim),
                 d_hessFD[ip](0, dim), d_hess[ip](1, dim), d_hessFD[ip](1, dim), d_hess[ip](2, dim),
                 d_hessFD[ip](2, dim));
 #endif
   }

◆ evaluateDerivativesFD()

bool evaluateDerivativesFD	(	const real_type	r_12,
		const real_type	r_1I,
		const real_type	r_2I,
		std::vector< double > &	d_vals,
		std::vector< TinyVector< real_type, 3 >> &	d_grads,
		std::vector< Tensor< real_type, 3 >> &	d_hess
	)

inline

Definition at line 703 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::evaluate(), PolynomialFunctor3D::Parameters, and PolynomialFunctor3D::reset_gamma().

Referenced by PolynomialFunctor3D::evaluateDerivatives().

   {
     const real_type eps = 1.0e-6;
     assert(d_vals.size() == Parameters.size());
     assert(d_grads.size() == Parameters.size());
     assert(d_hess.size() == Parameters.size());
     for (int ip = 0; ip < Parameters.size(); ip++)
     {
       real_type v_plus, v_minus;
       TinyVector<real_type, 3> g_plus, g_minus;
       Tensor<real_type, 3> h_plus, h_minus;
       real_type save_p = Parameters[ip];
       Parameters[ip]   = save_p + eps;
       reset_gamma();
       v_plus         = evaluate(r_12, r_1I, r_2I, g_plus, h_plus);
       Parameters[ip] = save_p - eps;
       reset_gamma();
       v_minus        = evaluate(r_12, r_1I, r_2I, g_minus, h_minus);
       Parameters[ip] = save_p;
       reset_gamma();
       real_type dp_inv = 0.5 / eps;
       d_vals[ip]       = dp_inv * (v_plus - v_minus);
       d_grads[ip]      = dp_inv * (g_plus - g_minus);
       d_hess[ip]       = dp_inv * (h_plus - h_minus);
     }
     return true;
   }

◆ evaluateV()

real_type evaluateV	(	int	Nptcl,
		const real_type *restrict	r_12_array,
		const real_type *restrict	r_1I_array,
		const real_type *restrict	r_2I_array
	)		const

inline

Definition at line 352 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::C, BLAS::cone, OptimizableFunctorBase::cutoff_radius, BLAS::czero, PolynomialFunctor3D::gamma, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, and PolynomialFunctor3D::N_eI.

   {
     constexpr real_type czero(0);
     constexpr real_type cone(1);
     constexpr real_type chalf(0.5);
 
     const real_type L = chalf * cutoff_radius;
     real_type val_tot = czero;
 
 #pragma omp simd aligned(r_12_array, r_1I_array, r_2I_array : QMC_SIMD_ALIGNMENT) reduction(+ : val_tot)
     for (int ptcl = 0; ptcl < Nptcl; ptcl++)
     {
       const real_type r_12 = r_12_array[ptcl];
       const real_type r_1I = r_1I_array[ptcl];
       const real_type r_2I = r_2I_array[ptcl];
       real_type val        = czero;
       real_type r2l(cone);
       for (int l = 0; l <= N_eI; l++)
       {
         real_type r2m(r2l);
         for (int m = 0; m <= N_eI; m++)
         {
           real_type r2n(r2m);
           for (int n = 0; n <= N_ee; n++)
           {
             val += gamma(l, m, n) * r2n;
             r2n *= r_12;
           }
           r2m *= r_2I;
         }
         r2l *= r_1I;
       }
       const real_type both_minus_L = (r_2I - L) * (r_1I - L);
       for (int i = 0; i < C; i++)
         val *= both_minus_L;
       val_tot += val;
     }
 
     return val_tot;
   }

◆ evaluateVGL()

void evaluateVGL	(	int	Nptcl,
		const real_type *restrict	r_12_array,
		const real_type *restrict	r_1I_array,
		const real_type *restrict	r_2I_array,
		real_type *restrict	val_array,
		real_type *restrict	grad0_array,
		real_type *restrict	grad1_array,
		real_type *restrict	grad2_array,
		real_type *restrict	hess00_array,
		real_type *restrict	hess11_array,
		real_type *restrict	hess22_array,
		real_type *restrict	hess01_array,
		real_type *restrict	hess02_array
	)		const

inline

Definition at line 481 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::C, BLAS::cone, OptimizableFunctorBase::cutoff_radius, BLAS::czero, PolynomialFunctor3D::gamma, qmcplusplus::Units::distance::m, qmcplusplus::n, PolynomialFunctor3D::N_ee, and PolynomialFunctor3D::N_eI.

   {
     constexpr real_type czero(0);
     constexpr real_type cone(1);
     constexpr real_type chalf(0.5);
     constexpr real_type ctwo(2);
 
     const real_type L = chalf * cutoff_radius;
 #pragma omp simd aligned(r_12_array, r_1I_array, r_2I_array, val_array, grad0_array, grad1_array, grad2_array, \
                          hess00_array, hess11_array, hess22_array, hess01_array, hess02_array                  \
                          : QMC_SIMD_ALIGNMENT)
     for (int ptcl = 0; ptcl < Nptcl; ptcl++)
     {
       const real_type r_12 = r_12_array[ptcl];
       const real_type r_1I = r_1I_array[ptcl];
       const real_type r_2I = r_2I_array[ptcl];
 
       real_type val(czero);
       real_type grad0(czero);
       real_type grad1(czero);
       real_type grad2(czero);
       real_type hess00(czero);
       real_type hess11(czero);
       real_type hess22(czero);
       real_type hess01(czero);
       real_type hess02(czero);
 
       real_type r2l(cone), r2l_1(czero), r2l_2(czero), lf(czero);
       for (int l = 0; l <= N_eI; l++)
       {
         real_type r2m(cone), r2m_1(czero), r2m_2(czero), mf(czero);
         for (int m = 0; m <= N_eI; m++)
         {
           real_type r2n(cone), r2n_1(czero), r2n_2(czero), nf(czero);
           for (int n = 0; n <= N_ee; n++)
           {
             const real_type g    = gamma(l, m, n);
             const real_type g00x = g * r2l * r2m;
             const real_type g10x = g * r2l_1 * r2m;
             const real_type g01x = g * r2l * r2m_1;
             const real_type gxx0 = g * r2n;
 
             val += g00x * r2n;
             grad0 += g00x * r2n_1;
             grad1 += g10x * r2n;
             grad2 += g01x * r2n;
             hess00 += g00x * r2n_2;
             hess01 += g10x * r2n_1;
             hess02 += g01x * r2n_1;
             hess11 += gxx0 * r2l_2 * r2m;
             hess22 += gxx0 * r2l * r2m_2;
             nf += cone;
             r2n_2 = r2n_1 * nf;
             r2n_1 = r2n * nf;
             r2n *= r_12;
           }
           mf += cone;
           r2m_2 = r2m_1 * mf;
           r2m_1 = r2m * mf;
           r2m *= r_2I;
         }
         lf += cone;
         r2l_2 = r2l_1 * lf;
         r2l_1 = r2l * lf;
         r2l *= r_1I;
       }
 
       const real_type r_2I_minus_L = r_2I - L;
       const real_type r_1I_minus_L = r_1I - L;
       const real_type both_minus_L = r_2I_minus_L * r_1I_minus_L;
       for (int i = 0; i < C; i++)
       {
         hess00 = both_minus_L * hess00;
         hess01 = both_minus_L * hess01 + r_2I_minus_L * grad0;
         hess02 = both_minus_L * hess02 + r_1I_minus_L * grad0;
         hess11 = both_minus_L * hess11 + ctwo * r_2I_minus_L * grad1;
         hess22 = both_minus_L * hess22 + ctwo * r_1I_minus_L * grad2;
         grad0  = both_minus_L * grad0;
         grad1  = both_minus_L * grad1 + r_2I_minus_L * val;
         grad2  = both_minus_L * grad2 + r_1I_minus_L * val;
         val *= both_minus_L;
       }
 
       val_array[ptcl]    = val;
       grad0_array[ptcl]  = grad0 / r_12;
       grad1_array[ptcl]  = grad1 / r_1I;
       grad2_array[ptcl]  = grad2 / r_2I;
       hess00_array[ptcl] = hess00;
       hess11_array[ptcl] = hess11;
       hess22_array[ptcl] = hess22;
       hess01_array[ptcl] = hess01 / (r_12 * r_1I);
       hess02_array[ptcl] = hess02 / (r_12 * r_2I);
     }
   }

◆ f()

real_type f ( real_type r )

inlineoverride

Definition at line 968 of file PolynomialFunctor3D.h.

968 { return 0.0; }

◆ getNumParameters()

int getNumParameters ( )

inline

Definition at line 1105 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::Parameters.

1105 { return Parameters.size(); }

qmcplusplus::PolynomialFunctor3D::Parameters

std::vector< real_type > Parameters

Definition: PolynomialFunctor3D.h:46

◆ makeClone()

OptimizableFunctorBase* makeClone ( ) const

inlineoverridevirtual

create a clone of this object

Implements OptimizableFunctorBase.

Definition at line 69 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::PolynomialFunctor3D().

69 { return new PolynomialFunctor3D(*this); }

qmcplusplus::PolynomialFunctor3D::PolynomialFunctor3D

PolynomialFunctor3D(const std::string &my_name, real_type ee_cusp=0.0, real_type eI_cusp=0.0)

constructor

Definition: PolynomialFunctor3D.h:58

◆ print()

void print ( std::ostream & os )

inline

Definition at line 1089 of file PolynomialFunctor3D.h.

   {
     /* no longer correct. Ye Luo
     int n=100;
     real_type d=cutoff_radius/100.,r=0;
     real_type u,du,d2du;
     for(int i=0; i<n; ++i)
     {
       u=evaluate(r,du,d2du);
       os << std::setw(22) << r << std::setw(22) << u << std::setw(22) << du
          << std::setw(22) << d2du << std::endl;
       r+=d;
     }
     */
   }

◆ put()

bool put ( xmlNodePtr cur )

inlineoverridevirtual

process xmlnode and registers variables to optimize

Parameters

cur	xmlNode for a functor

Implements OptimizableFunctorBase.

Definition at line 971 of file PolynomialFunctor3D.h.

   {
     ReportEngine PRE("PolynomialFunctor3D", "put(xmlNodePtr)");
     // //CuspValue = -1.0e10;
     // NumParams_eI = NumParams_ee = 0;
     cutoff_radius = 0.0;
     OhmmsAttributeSet rAttrib;
     rAttrib.add(N_ee, "esize");
     rAttrib.add(N_eI, "isize");
     rAttrib.add(cutoff_radius, "rcut");
     rAttrib.put(cur);
     if (N_eI == 0)
       PRE.error("You must specify a positive number for \"isize\"", true);
     if (N_ee == 0)
       PRE.error("You must specify a positive number for \"esize\"", true);
     app_summary() << "     Ion: " << iSpecies << "   electron-electron: " << eSpecies1 << " - " << eSpecies2
                   << std::endl;
     app_summary() << "      Number of parameters for e-e: " << N_ee << ", for e-I: " << N_eI << std::endl;
     app_summary() << "      Cutoff radius: " << cutoff_radius << std::endl;
     app_summary() << std::endl;
     resize(N_eI, N_ee);
     // Now read coefficents
     xmlNodePtr xmlCoefs = cur->xmlChildrenNode;
     while (xmlCoefs != NULL)
     {
       std::string cname((const char*)xmlCoefs->name);
       if (cname == "coefficients")
       {
         std::string type("0"), id("0"), opt("yes");
         OhmmsAttributeSet cAttrib;
         cAttrib.add(id, "id");
         cAttrib.add(type, "type");
         cAttrib.add(opt, "optimize");
         cAttrib.put(xmlCoefs);
         notOpt = (opt == "no");
         if (type != "Array")
         {
           PRE.error("Unknown correlation type " + type + " in PolynomialFunctor3D." + "Resetting to \"Array\"");
           xmlNewProp(xmlCoefs, (const xmlChar*)"type", (const xmlChar*)"Array");
         }
         std::vector<real_type> params;
         putContent(params, xmlCoefs);
         if (params.size() == Parameters.size())
           Parameters = params;
         else
         {
           app_log() << "Expected " << Parameters.size() << " parameters,"
                     << " but found " << params.size() << " in PolynomialFunctor3D.\n";
           if (params.size() != 0)
             abort(); //you think you know what they should be but don't.
         }
         // Setup parameter names
         for (int i = 0; i < Parameters.size(); i++)
         {
           std::stringstream sstr;
           sstr << id << "_" << i;
           ;
           if (!notOpt)
             myVars.insert(sstr.str(), Parameters[i], optimize::LOGLINEAR_P, true);
         }
         // for (int i=0; i< N_ee; i++)
         //   for (int j=0; j < N_eI; j++)
         //     for (int k=0; k<=j; k++) {
         //  std::stringstream sstr;
         //  sstr << id << "_" << i << "_" << j << "_" << k;
         //  myVars.insert(sstr.str(),Parameters[index],true);
         //  ParamArray(i,j,k) = ParamArray(i,k,j) = Parameters[index];
         //  index++;
         //     }
         if (!notOpt)
         {
           int left_pad_spaces = 6;
           myVars.print(app_log(), left_pad_spaces, true);
           app_log() << std::endl;
         }
       }
       xmlCoefs = xmlCoefs->next;
     }
     reset_gamma();
     return true;
   }

◆ reset()

void reset ( )

inlineoverridevirtual

reset function

Implements OptimizableFunctorBase.

Definition at line 219 of file PolynomialFunctor3D.h.

References PolynomialFunctor3D::N_ee, PolynomialFunctor3D::N_eI, PolynomialFunctor3D::reset_gamma(), and PolynomialFunctor3D::resize().

   {
     resize(N_eI, N_ee);
     reset_gamma();
   }

◆ reset_gamma()

void reset_gamma ( )

inline

Definition at line 225 of file PolynomialFunctor3D.h.

Referenced by PolynomialFunctor3D::evaluateDerivativesFD(), PolynomialFunctor3D::put(), PolynomialFunctor3D::reset(), and PolynomialFunctor3D::resetParametersExclusive().

   {
     // fprintf (stderr, "Parameters:\n");
     // for (int i=0; i<Parameters.size(); i++)
     //  fprintf (stderr, " %16.10e\n", Parameters[i]);
     const double L = 0.5 * cutoff_radius;
     std::fill(GammaVec.begin(), GammaVec.end(), 0.0);
     // First, set all independent variables
     int var = 0;
     for (int i = 0; i < NumGamma; i++)
       if (IndepVar[i])
         GammaVec[i] = scale * Parameters[var++];
     assert(var == Parameters.size());
     // Now, set dependent variables
     var = 0;
     //      std::cerr << "NumConstraints = " << NumConstraints << std::endl;
     for (int i = 0; i < NumGamma; i++)
       if (!IndepVar[i])
       {
         // fprintf (stderr, "constraintMatrix(%d,%d) = %1.10f\n",
         //     var, i, ConstraintMatrix(var,i));
         assert(std::abs(ConstraintMatrix(var, i) - 1.0) < 1.0e-6);
         for (int j = 0; j < NumGamma; j++)
           if (i != j)
             GammaVec[i] -= ConstraintMatrix(var, j) * GammaVec[j];
         var++;
       }
     int num = 0;
     for (int m = 0; m <= N_eI; m++)
       for (int l = m; l <= N_eI; l++)
         for (int n = 0; n <= N_ee; n++)
           //      gamma(m,l,n) = gamma(l,m,n) = unpermuted[num++];
           gamma(m, l, n) = gamma(l, m, n) = GammaVec[num++];
     // Now check that constraints have been satisfied
     // e-e constraints
     for (int k = 0; k <= 2 * N_eI; k++)
     {
       real_type sum = 0.0;
       for (int m = 0; m <= k; m++)
       {
         int l = k - m;
         if (l <= N_eI && m <= N_eI)
         {
           int i = index(l, m, 1);
           if (l > m)
             sum += 2.0 * GammaVec[i];
           else if (l == m)
             sum += GammaVec[i];
         }
       }
       if (std::abs(sum) > 1.0e-9)
         std::cerr << "error in k = " << k << "  sum = " << sum << std::endl;
     }
     for (int k = 0; k <= 2 * N_eI; k++)
     {
       real_type sum = 0.0;
       for (int l = 0; l <= k; l++)
       {
         int m = k - l;
         if (m <= N_eI && l <= N_eI)
         {
           // fprintf (stderr, "k = %d gamma(%d, %d, 1) = %1.8f\n", k, l, m,
           //       gamma(l,m,1));
           sum += gamma(l, m, 1);
         }
       }
       if (std::abs(sum) > 1.0e-6)
       {
         app_error() << "e-e constraint not satisfied in PolynomialFunctor3D:  k=" << k << "  sum=" << sum << std::endl;
         abort();
       }
     }
     // e-I constraints
     for (int k = 0; k <= N_eI + N_ee; k++)
     {
       real_type sum = 0.0;
       for (int m = 0; m <= k; m++)
       {
         int n = k - m;
         if (m <= N_eI && n <= N_ee)
         {
           sum += (real_type)C * gamma(0, m, n) - L * gamma(1, m, n);
           // fprintf (stderr,
           //       "k = %d gamma(0,%d,%d) = %1.8f  gamma(1,%d,%d)=%1.8f\n",
           //       k, m, n, gamma(0,m,n), m, n, gamma(1,m,n));
         }
       }
       if (std::abs(sum) > 1.0e-6)
       {
         app_error() << "e-I constraint not satisfied in PolynomialFunctor3D:  k=" << k << "  sum=" << sum << std::endl;
         abort();
       }
     }
   }

◆ resetParametersExclusive()

void resetParametersExclusive ( const opt_variables_type & active )

inlineoverridevirtual

reset the parameters during optimizations.

Exclusive, see checkInVariablesExclusive

Implements OptimizableObject.

Definition at line 1053 of file PolynomialFunctor3D.h.

References OptimizableFunctorBase::myVars, PolynomialFunctor3D::notOpt, PolynomialFunctor3D::Parameters, PolynomialFunctor3D::reset_gamma(), PolynomialFunctor3D::ResetCount, and VariableSet::where().

   {
     if (notOpt)
       return;
 
     for (int i = 0; i < Parameters.size(); ++i)
     {
       int loc = myVars.where(i);
       if (loc >= 0)
       {
         Parameters[i] = std::real(myVars[i] = active[loc]);
       }
     }
 
     if (ResetCount++ == 100)
       ResetCount = 0;
 
     reset_gamma();
   }

◆ resize()

void resize	(	int	neI,
		int	nee
	)

inline

Definition at line 71 of file PolynomialFunctor3D.h.

Referenced by PolynomialFunctor3D::put(), and PolynomialFunctor3D::reset().

   {
     N_eI           = neI;
     N_ee           = nee;
     const double L = 0.5 * cutoff_radius;
     gamma.resize(N_eI + 1, N_eI + 1, N_ee + 1);
     index.resize(N_eI + 1, N_eI + 1, N_ee + 1);
     NumGamma       = ((N_eI + 1) * (N_eI + 2) / 2 * (N_ee + 1));
     NumConstraints = (2 * N_eI + 1) + (N_eI + N_ee + 1);
     int numParams  = NumGamma - NumConstraints;
     Parameters.resize(numParams);
     d_valsFD.resize(numParams);
     d_gradsFD.resize(numParams);
     d_hessFD.resize(numParams);
     GammaVec.resize(NumGamma);
     dval_Vec.resize(NumGamma);
     dgrad_Vec.resize(NumGamma);
     dhess_Vec.resize(NumGamma);
     ConstraintMatrix.resize(NumConstraints, NumGamma);
     // Assign indices
     int num = 0;
     for (int m = 0; m <= N_eI; m++)
       for (int l = m; l <= N_eI; l++)
         for (int n = 0; n <= N_ee; n++)
           index(l, m, n) = index(m, l, n) = num++;
     assert(num == NumGamma);
     //       std::cerr << "NumGamma = " << NumGamma << std::endl;
     // Fill up contraint matrix
     // For 3 constraints and 2 parameters, we would have
     // |A00 A01 A02 A03 A04|  |g0|   |0 |
     // |A11 A11 A12 A13 A14|  |g1|   |0 |
     // |A22 A21 A22 A23 A24|  |g2| = |0 |
     // | 0   0   0   1   0 |  |g3|   |p0|
     // | 0   0   0   0   1 |  |g4|   |p1|
     ConstraintMatrix = 0.0;
     // std::cerr << "ConstraintMatrix.size = " << ConstraintMatrix.size(0)
     //     << " by " << ConstraintMatrix.size(1) << std::endl;
     // std::cerr << "index.size() = (" << index.size(0) << ", "
     //     << index.size(1) << ", " << index.size(2) << ").\n";
     int k;
     // e-e no-cusp constraint
     for (k = 0; k <= 2 * N_eI; k++)
     {
       for (int m = 0; m <= k; m++)
       {
         int l = k - m;
         if (l <= N_eI && m <= N_eI)
         {
           int i = index(l, m, 1);
           if (l > m)
             ConstraintMatrix(k, i) = 2.0;
           else if (l == m)
             ConstraintMatrix(k, i) = 1.0;
         }
       }
     }
     // e-I no-cusp constraint
     for (int kp = 0; kp <= N_eI + N_ee; kp++)
     {
       if (kp <= N_ee)
       {
         ConstraintMatrix(k + kp, index(0, 0, kp)) = (real_type)C;
         ConstraintMatrix(k + kp, index(0, 1, kp)) = -L;
       }
       for (int l = 1; l <= kp; l++)
       {
         int n = kp - l;
         if (n >= 0 && n <= N_ee && l <= N_eI)
         {
           ConstraintMatrix(k + kp, index(l, 0, n)) = (real_type)C;
           ConstraintMatrix(k + kp, index(l, 1, n)) = -L;
         }
       }
     }
     //    {
     //       fprintf (stderr, "Constraint matrix:\n");
     //       for (int i=0; i<NumConstraints; i++) {
     //  for (int j=0; j<NumGamma; j++)
     //    fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j));
     //  fprintf(stderr, "\n");
     //       }
     //     }
     // Now, row-reduce constraint matrix
     GammaPerm.resize(NumGamma);
     IndepVar.resize(NumGamma, false);
     // Set identity permutation
     for (int i = 0; i < NumGamma; i++)
       GammaPerm[i] = i;
     int col = -1;
     for (int row = 0; row < NumConstraints; row++)
     {
       int max_loc;
       real_type max_abs;
       do
       {
         col++;
         max_loc = row;
         max_abs = std::abs(ConstraintMatrix(row, col));
         for (int ri = row + 1; ri < NumConstraints; ri++)
         {
           real_type abs_val = std::abs(ConstraintMatrix(ri, col));
           if (abs_val > max_abs)
           {
             max_loc = ri;
             max_abs = abs_val;
           }
         }
         if (max_abs < 1.0e-6)
           IndepVar[col] = true;
       } while (max_abs < 1.0e-6);
       ConstraintMatrix.swap_rows(row, max_loc);
       real_type lead_inv = 1.0 / ConstraintMatrix(row, col);
       for (int c = 0; c < NumGamma; c++)
         ConstraintMatrix(row, c) *= lead_inv;
       // Now, eliminate column entries
       for (int ri = 0; ri < NumConstraints; ri++)
       {
         if (ri != row)
         {
           real_type val = ConstraintMatrix(ri, col);
           for (int c = 0; c < NumGamma; c++)
             ConstraintMatrix(ri, c) -= val * ConstraintMatrix(row, c);
         }
       }
     }
     for (int c = col + 1; c < NumGamma; c++)
       IndepVar[c] = true;
     //       fprintf (stderr, "Reduced Constraint matrix:\n");
     //       for (int i=0; i<NumConstraints; i++) {
     //  for (int j=0; j<NumGamma; j++)
     //    fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j));
     //  fprintf(stderr, "\n");
     //       }
     //       fprintf (stderr, "Independent vars = \n");
     //       for (int i=0; i<NumGamma; i++)
     //  if (IndepVar[i])
     //    fprintf (stderr, "%d ", i);
     //       fprintf (stderr, "\n");
     // fprintf (stderr, "Inverse matrix:\n");
     // // Now, invert constraint matrix
     // Invert(ConstraintMatrix.data(), NumGamma, NumGamma);
     // for (int i=0; i<NumGamma; i++) {
     //  for (int j=0; j<NumGamma; j++)
     //    fprintf (stderr, "%5.2f ", ConstraintMatrix(i,j));
     //  fprintf(stderr, "\n");
     // }
   }