One and two-body jastrows give Hessians wrt electron coordinates

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@6265 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Raymond Clay 2014-03-13 18:22:09 +00:00
parent cc7e97f6d9
commit 69a8b77541
2 changed files with 66 additions and 1 deletions

View File

@ -89,6 +89,7 @@ public:
typedef FT FuncType;
///constructor
OneBodyJastrowOrbital(const ParticleSet& centers, ParticleSet& els)
: CenterRef(centers), FirstAddressOfdU(0), LastAddressOfdU(0)
@ -223,6 +224,34 @@ public:
{
return std::exp(evaluateLog(P,G,L));
}
void evaluateHessian(ParticleSet& P, HessVector_t& grad_grad_psi)
{
LogValue=0.0;
U=0.0;
const DistanceTableData* d_table=P.DistTables[myTableIndex];
RealType dudr, d2udr2;
Tensor<RealType,OHMMS_DIM> ident;
grad_grad_psi=0.0;
ident.diagonal(1.0);
for (int i=0; i<d_table->size(SourceIndex); i++)
{
FT* func=Fs[i];
if (func == 0)
continue;
for (int nn=d_table->M[i]; nn<d_table->M[i+1]; nn++)
{
int j = d_table->J[nn];
RealType rinv=d_table->rinv(nn);
RealType uij= func->evaluate(d_table->r(nn), dudr, d2udr2);
grad_grad_psi[j]-= rinv*rinv*outerProduct(d_table->dr(nn),d_table->dr(nn))*(d2udr2-dudr*rinv) + ident*dudr*rinv;
}
}
}
// ValueType evaluate(ParticleSet& P, ParticleSet::
/** evaluate the ratio \f$exp(U(iat)-U_0(iat))\f$
* @param P active particle set

View File

@ -276,7 +276,43 @@ public:
{
return std::exp(evaluateLog(P,G,L));
}
void evaluateHessian(ParticleSet& P, HessVector_t& grad_grad_psi)
{
LogValue=0.0;
const DistanceTableData* d_table=P.DistTables[0];
RealType dudr, d2udr2;
PosType gr;
Tensor<RealType,OHMMS_DIM> ident;
grad_grad_psi=0.0;
ident.diagonal(1.0);
for(int i=0; i<d_table->size(SourceIndex); i++)
{
for(int nn=d_table->M[i]; nn<d_table->M[i+1]; nn++)
{
int j = d_table->J[nn];
//LogValue -= F[d_table->PairID[nn]]->evaluate(d_table->r(nn), dudr, d2udr2);
RealType uij = F[d_table->PairID[nn]]->evaluate(d_table->r(nn), dudr, d2udr2);
LogValue -= uij;
// U[i*N+j]=uij;
// U[j*N+i]=uij; //save for the ratio
//multiply 1/r
// dudr *= d_table->rinv(nn);
// gr = dudr*d_table->dr(nn);
//(d^2 u \over dr^2) + (2.0\over r)(du\over\dr)
RealType rinv = d_table->rinv(nn);
Tensor<RealType, OHMMS_DIM> hess = rinv*rinv*outerProduct(d_table->dr(nn),d_table->dr(nn))*(d2udr2-dudr*rinv) + ident*dudr*rinv;
grad_grad_psi[i] -= hess;
grad_grad_psi[j] -= hess;
}
}
}
ValueType ratio(ParticleSet& P, int iat)
{
const DistanceTableData* d_table=P.DistTables[0];