mirror of https://github.com/QMCPACK/qmcpack.git
complex works fully. real has hderiv issues with 3x3x3 tiling.
git-svn-id: https://subversion.assembla.com/svn/qmcdev/branches/OptBF@5456 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Jeremy McMinis
parent
63a43fefc4
commit
0140365fbd
|
|
@ -250,13 +250,30 @@ outerProduct(const TinyVector<T1,D> &lhs, const TinyVector<T2,D> &rhs)
|
|||
return OuterProduct< TinyVector<T1,D> , TinyVector<T2,D> > :: apply(lhs,rhs);
|
||||
}
|
||||
|
||||
template < class T1, unsigned D >
|
||||
inline TinyVector<Tensor<T1,D>,D>
|
||||
outerdot(const TinyVector<T1,D> &lhs, const TinyVector<T1,D> &mhs, const TinyVector<T1,D> &rhs)
|
||||
{
|
||||
TinyVector<Tensor<T1,D>,D> ret;
|
||||
Tensor<T1,D> tmp=OuterProduct< TinyVector<T1,D> , TinyVector<T1,D> > :: apply(lhs,mhs);
|
||||
for(unsigned i(0);i<D;i++) ret[i]=rhs[i]*tmp;
|
||||
return ret;
|
||||
}
|
||||
|
||||
template < class T1, class T2, class T3, unsigned D >
|
||||
inline TinyVector<Tensor<typename BinaryReturn<T1,T2,OpMultiply>::Type_t,D>,D>
|
||||
outerProduct(const TinyVector<T1,D> &lhs, const TinyVector<T2,D> &mhs, const TinyVector<T3,D> &rhs)
|
||||
symouterdot(const TinyVector<T1,D> &lhs, const TinyVector<T2,D> &mhs, const TinyVector<T3,D> &rhs)
|
||||
{
|
||||
TinyVector<Tensor<typename BinaryReturn<T1,T2,OpMultiply>::Type_t,D>,D> ret;
|
||||
Tensor<typename BinaryReturn<T1,T2,OpMultiply>::Type_t,D> tmp=OuterProduct< TinyVector<T1,D> , TinyVector<T2,D> > :: apply(lhs,mhs);
|
||||
for(unsigned i(0);i<D;i++) ret[i]=rhs[i]*tmp;
|
||||
|
||||
tmp=OuterProduct< TinyVector<T2,D> , TinyVector<T3,D> > :: apply(mhs,rhs);
|
||||
for(unsigned i(0);i<D;i++) ret[i]+=lhs[i]*tmp;
|
||||
|
||||
tmp=OuterProduct< TinyVector<T1,D> , TinyVector<T3,D> > :: apply(lhs,rhs);
|
||||
for(unsigned i(0);i<D;i++) ret[i]+=mhs[i]*tmp;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1758,152 +1758,6 @@ namespace qmcplusplus {
|
|||
RealGGGMatrix_t& grad_grad_grad_logdet)
|
||||
{
|
||||
// APP_ABORT(" EinsplineSetExtended<StorageType>::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n");
|
||||
|
||||
complex<double> eye(0.0,1.0);
|
||||
VGLMatTimer.start();
|
||||
for (int iat=first,i=0; iat<last; iat++,i++) {
|
||||
PosType r (P.R[iat]);
|
||||
|
||||
// Do core states first
|
||||
int icore = NumValenceOrbs;
|
||||
for (int tin=0; tin<MuffinTins.size(); tin++) {
|
||||
APP_ABORT("MuffinTins not implemented with Hessian evaluation.\n");
|
||||
}
|
||||
|
||||
// Check if we are in the muffin tin; if so, evaluate
|
||||
bool inTin = false, need2blend = false;
|
||||
PosType disp;
|
||||
double b, db, d2b;
|
||||
for (int tin=0; tin<MuffinTins.size(); tin++) {
|
||||
APP_ABORT("MuffinTins not implemented with Hessian evaluation.\n");
|
||||
}
|
||||
|
||||
|
||||
bool inAtom = false;
|
||||
for (int jat=0; jat<AtomicOrbitals.size(); jat++) {
|
||||
inAtom = AtomicOrbitals[jat].evaluate
|
||||
(r, StorageValueVector, StorageGradVector, StorageHessVector);
|
||||
if (inAtom) break;
|
||||
}
|
||||
|
||||
StorageValueVector_t &valVec =
|
||||
need2blend ? BlendValueVector : StorageValueVector;
|
||||
StorageGradVector_t &gradVec =
|
||||
need2blend ? BlendGradVector : StorageGradVector;
|
||||
StorageHessVector_t &hessVec =
|
||||
need2blend ? BlendHessVector : StorageHessVector;
|
||||
StorageGradHessVector_t &gradhessVec = StorageGradHessVector;
|
||||
|
||||
Tensor<complex<double>,OHMMS_DIM> tmphs;
|
||||
TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
// Otherwise, evaluate the B-splines
|
||||
if (!inTin || need2blend) {
|
||||
if (!inAtom) {
|
||||
PosType ru(PrimLattice.toUnit(P.R[iat]));
|
||||
for (int i=0; i<OHMMS_DIM; i++)
|
||||
ru[i] -= std::floor (ru[i]);
|
||||
EinsplineTimer.start();
|
||||
// EinsplineMultiEval (MultiSpline, ru, valVec, gradVec, StorageHessVector);
|
||||
EinsplineMultiEval (MultiSpline, ru, valVec, gradVec, StorageHessVector,StorageGradHessVector);
|
||||
EinsplineTimer.stop();
|
||||
Tensor<double,OHMMS_DIM> TPG(transpose(PrimLattice.G));
|
||||
Tensor<double,OHMMS_DIM> PG2;
|
||||
// for (int n=0; n<OHMMS_DIM; n++) for (int m=0; m<OHMMS_DIM; m++) PG2(n,m)=TPG(n,m)*TPG(n,m);
|
||||
PG2=dot(TPG,TPG);
|
||||
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
gradVec[j] = dot (PrimLattice.G, gradVec[j]);
|
||||
// FIX FIX FIX: store transpose(PrimLattice.G)
|
||||
tmphs = dot(TPG,StorageHessVector[j]);
|
||||
hessVec[j] = dot(tmphs,PrimLattice.G);
|
||||
//Is this right?
|
||||
tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
gradhessVec[j]=dot(tmpghs,PrimLattice.G);
|
||||
// gradhessVec[j]=dot(StorageGradHessVector[j],PrimLattice.G);
|
||||
}
|
||||
|
||||
// Add e^-ikr phase to B-spline orbitals
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
complex<double> u = valVec[j];
|
||||
TinyVector<complex<double>,OHMMS_DIM> gradu = gradVec[j];
|
||||
tmphs = hessVec[j];
|
||||
tmpghs = gradhessVec[j];
|
||||
PosType k = kPoints[j];
|
||||
TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
for (int n=0; n<OHMMS_DIM; n++) ck[n] = k[n];
|
||||
double s,c;
|
||||
double phase = -dot(r, k);
|
||||
sincos (phase, &s, &c);
|
||||
complex<double> e_mikr (c,s);
|
||||
valVec[j] = e_mikr*u;
|
||||
gradVec[j] = e_mikr*(-eye*u*ck + gradu);
|
||||
hessVec[j] = e_mikr*(tmphs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck));
|
||||
//Is this right?
|
||||
for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
|
||||
gradhessVec[j]=e_mikr*(tmpghs +eye*u*outerProduct(ck,ck,ck)
|
||||
-outerProduct(gradu,ck,ck) -outerProduct(ck,gradu,ck) -outerProduct(ck,ck,gradu) -hvdot );
|
||||
}
|
||||
}
|
||||
}
|
||||
// Finally, copy into output vectors
|
||||
int psiIndex = 0;
|
||||
int N = StorageValueVector.size();
|
||||
if (need2blend) {
|
||||
|
||||
APP_ABORT("need2blend not implemented with Hessian evaluation.\n");
|
||||
}
|
||||
else { // No blending needed
|
||||
for (int j=0; j<N; j++) {
|
||||
// complex<double> psi_val;
|
||||
// TinyVector<complex<double>, OHMMS_DIM> psi_grad;
|
||||
// psi_val = StorageValueVector[j];
|
||||
// psi_grad = StorageGradVector[j];
|
||||
// tmphs = StorageHessVector[j];
|
||||
// tmpghs = StorageGradHessVector[j];
|
||||
//
|
||||
psi(i,psiIndex)=real(StorageValueVector[j]);
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
dpsi(i,psiIndex)(n) = real(StorageGradVector[j](n));
|
||||
for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
grad_grad_psi(i,psiIndex)(n) = real(StorageHessVector[j](n));
|
||||
////Is this right?
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
grad_grad_grad_logdet(i,psiIndex)[n](m) = real(StorageGradHessVector[j][n](m));
|
||||
psiIndex++;
|
||||
// // if (psiIndex >= dpsi.cols()) {
|
||||
// // cerr << "Error: out of bounds writing in EinsplineSet::evalate.\n"
|
||||
// // << "psiIndex = " << psiIndex << " dpsi.cols() = " << dpsi.cols() << endl;
|
||||
// // }
|
||||
if (MakeTwoCopies[j]) {
|
||||
psi(i,psiIndex)=imag(StorageValueVector[j]);
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
dpsi(i,psiIndex)(n) = imag(StorageGradVector[j](n));
|
||||
for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
grad_grad_psi(i,psiIndex)(n) = imag(StorageGradVector[j](n));
|
||||
////Is this right?
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
grad_grad_grad_logdet(i,psiIndex)[n](m) = imag(StorageGradHessVector[j][n](m));
|
||||
psiIndex++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
VGLMatTimer.stop();
|
||||
}
|
||||
|
||||
template<> void
|
||||
EinsplineSetExtended<double>::evaluate_notranspose(const ParticleSet& P, int first, int last,
|
||||
ComplexValueMatrix_t& psi, ComplexGradMatrix_t& dpsi,
|
||||
ComplexHessMatrix_t& grad_grad_psi,
|
||||
ComplexGGGMatrix_t& grad_grad_grad_logdet)
|
||||
{
|
||||
// APP_ABORT(" EinsplineSetExtended<StorageType>::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n");
|
||||
VGLMatTimer.start();
|
||||
for(int iat=first,i=0; iat<last; iat++,i++) {
|
||||
PosType r (P.R[iat]);
|
||||
|
|
@ -1911,52 +1765,317 @@ namespace qmcplusplus {
|
|||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
ru[n] -= std::floor (ru[n]);
|
||||
EinsplineTimer.start();
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
|
||||
StorageGradVector, StorageHessVector,StorageGradHessVector);
|
||||
EinsplineTimer.stop();
|
||||
Tensor<double,OHMMS_DIM> TPG(transpose(PrimLattice.G));
|
||||
Tensor<double,OHMMS_DIM> PG2;
|
||||
// for (int n=0; n<OHMMS_DIM; n++) for (int m=0; m<OHMMS_DIM; m++) PG2(n,m)=TPG(n,m)*TPG(n,m);
|
||||
PG2=dot(TPG,TPG);
|
||||
TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
complex<double> eye (0.0, 1.0);
|
||||
for (int j=0; j<OrbitalSetSize; j++) {
|
||||
complex<double> u;
|
||||
TinyVector<complex<double>, OHMMS_DIM> gradu;
|
||||
Tensor<complex<double>,OHMMS_DIM> hs,tmphs;
|
||||
Tensor<complex<double>,OHMMS_DIM> PG; PG=PrimLattice.G;
|
||||
Tensor<complex<double>,OHMMS_DIM> TPG; TPG=transpose(PrimLattice.G);
|
||||
Tensor<complex<double>,OHMMS_DIM> PG2;
|
||||
|
||||
u = StorageValueVector[j];
|
||||
gradu = dot(PrimLattice.G, StorageGradVector[j]);
|
||||
tmphs = dot(transpose(PrimLattice.G),StorageHessVector[j]);
|
||||
hs = dot(tmphs,PrimLattice.G);
|
||||
tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
tmpghs = dot(tmpghs,PrimLattice.G);
|
||||
PG2=dot(TPG,TPG);
|
||||
Tensor<complex<double>,OHMMS_DIM> hs,tmphs;
|
||||
TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
StorageGradVector[j] = dot (PG, StorageGradVector[j]);
|
||||
tmphs = dot(TPG,StorageHessVector[j]);
|
||||
StorageHessVector[j] = dot(tmphs,PG);
|
||||
tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
StorageGradHessVector[j]=dot(tmpghs,PG);
|
||||
}
|
||||
|
||||
complex<double> eye (0.0, 1.0);
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
complex<double> u(StorageValueVector[j]);
|
||||
TinyVector<complex<double>,OHMMS_DIM> gradu(StorageGradVector[j]);
|
||||
tmphs = StorageHessVector[j];
|
||||
tmpghs = StorageGradHessVector[j];
|
||||
|
||||
PosType k = kPoints[j];
|
||||
TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
ck[n] = k[n];
|
||||
for (int n=0; n<OHMMS_DIM; n++) ck[n] = k[n];
|
||||
double s,c;
|
||||
double phase = -dot(r, k);
|
||||
sincos (phase, &s, &c);
|
||||
complex<double> e_mikr (c,s);
|
||||
convert(e_mikr * u, psi(i,j));
|
||||
//convert(e_mikr * u, psi(j,i));
|
||||
convert(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
|
||||
//convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
|
||||
convert(e_mikr*(hs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck)),grad_grad_psi(i,j));
|
||||
for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
|
||||
convert(e_mikr*(tmpghs +eye*u*outerProduct(ck,ck,ck)
|
||||
-outerProduct(gradu,ck,ck) -outerProduct(ck,gradu,ck) -outerProduct(ck,ck,gradu) -hvdot ), grad_grad_grad_logdet(i,j));
|
||||
StorageValueVector[j] = e_mikr*u;
|
||||
StorageGradVector[j] = e_mikr*(-eye*u*ck + gradu);
|
||||
StorageHessVector[j] = e_mikr*(tmphs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck));
|
||||
//Is this right?
|
||||
for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
|
||||
StorageGradHessVector[j]=e_mikr*(tmpghs +eye*u*outerdot(ck,ck,ck) -symouterdot(gradu,ck,ck) -hvdot );
|
||||
}
|
||||
int psiIndex(0);
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
cerr<<psiIndex<<" "<<i<<" "<<j<<endl;
|
||||
psi(i,psiIndex)=real(StorageValueVector[j]);
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
dpsi(i,psiIndex)[n] = real(StorageGradVector[j][n]);
|
||||
for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
grad_grad_psi(i,psiIndex)(n) = real(StorageHessVector[j](n));
|
||||
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
grad_grad_grad_logdet(i,psiIndex)[n](m) = real(StorageGradHessVector[j][n](m));
|
||||
psiIndex++;
|
||||
|
||||
if (MakeTwoCopies[j]) {
|
||||
psi(i,psiIndex)=imag(StorageValueVector[j]);
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
dpsi(i,psiIndex)[n] = imag(StorageGradVector[j][n]);
|
||||
for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
grad_grad_psi(i,psiIndex)(n) = imag(StorageGradVector[j](n));
|
||||
|
||||
for (int n=0; n<OHMMS_DIM; n++)
|
||||
for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
grad_grad_grad_logdet(i,psiIndex)[n](m) = imag(StorageGradHessVector[j][n](m));
|
||||
psiIndex++;
|
||||
}
|
||||
}
|
||||
}
|
||||
VGLMatTimer.stop();
|
||||
}
|
||||
|
||||
|
||||
|
||||
// template<typename StorageType> void
|
||||
// EinsplineSetExtended<StorageType>::evaluate_notranspose(const ParticleSet& P, int first, int last,
|
||||
// RealValueMatrix_t& psi, RealGradMatrix_t& dpsi,
|
||||
// RealHessMatrix_t& grad_grad_psi,
|
||||
// RealGGGMatrix_t& grad_grad_grad_logdet)
|
||||
// {
|
||||
//// APP_ABORT(" EinsplineSetExtended<StorageType>::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n");
|
||||
//
|
||||
// complex<double> eye(0.0,1.0);
|
||||
// VGLMatTimer.start();
|
||||
// for (int iat=first,i=0; iat<last; iat++,i++) {
|
||||
// PosType r (P.R[iat]);
|
||||
//
|
||||
// // Do core states first
|
||||
// int icore = NumValenceOrbs;
|
||||
// for (int tin=0; tin<MuffinTins.size(); tin++) {
|
||||
// APP_ABORT("MuffinTins not implemented with Hessian evaluation.\n");
|
||||
// }
|
||||
//
|
||||
// // Check if we are in the muffin tin; if so, evaluate
|
||||
// bool inTin = false, need2blend = false;
|
||||
// PosType disp;
|
||||
// double b, db, d2b;
|
||||
// for (int tin=0; tin<MuffinTins.size(); tin++) {
|
||||
// APP_ABORT("MuffinTins not implemented with Hessian evaluation.\n");
|
||||
// }
|
||||
//
|
||||
//
|
||||
// bool inAtom = false;
|
||||
// for (int jat=0; jat<AtomicOrbitals.size(); jat++) {
|
||||
// inAtom = AtomicOrbitals[jat].evaluate
|
||||
// (r, StorageValueVector, StorageGradVector, StorageHessVector);
|
||||
// if (inAtom) break;
|
||||
// }
|
||||
//
|
||||
// StorageValueVector_t &valVec =
|
||||
// need2blend ? BlendValueVector : StorageValueVector;
|
||||
// StorageGradVector_t &gradVec =
|
||||
// need2blend ? BlendGradVector : StorageGradVector;
|
||||
// StorageHessVector_t &hessVec =
|
||||
// need2blend ? BlendHessVector : StorageHessVector;
|
||||
// StorageGradHessVector_t &gradhessVec = StorageGradHessVector;
|
||||
//
|
||||
// Tensor<complex<double>,OHMMS_DIM> tmphs;
|
||||
// TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
// // Otherwise, evaluate the B-splines
|
||||
// if (!inTin || need2blend) {
|
||||
// if (!inAtom) {
|
||||
// PosType ru(PrimLattice.toUnit(P.R[iat]));
|
||||
// for (int i=0; i<OHMMS_DIM; i++)
|
||||
// ru[i] -= std::floor (ru[i]);
|
||||
//
|
||||
//// int sign=0;
|
||||
//// for (int n=0; n<OHMMS_DIM; n++) {
|
||||
//// RealType img = std::floor(ru[n]);
|
||||
//// ru[n] -= img;
|
||||
//// sign += HalfG[n] * (int)img;
|
||||
//// }
|
||||
//
|
||||
//
|
||||
// EinsplineTimer.start();
|
||||
//// EinsplineMultiEval (MultiSpline, ru, valVec, gradVec, StorageHessVector);
|
||||
// EinsplineMultiEval (MultiSpline, ru, valVec, gradVec, StorageHessVector,StorageGradHessVector);
|
||||
// EinsplineTimer.stop();
|
||||
//
|
||||
//// if (sign & 1)
|
||||
//// for (int j=0; j<NumValenceOrbs; j++) {
|
||||
//// StorageValueVector[j] *= -1.0;
|
||||
//// StorageGradVector[j] *= -1.0;
|
||||
//// StorageHessVector[j] *= -1.0;
|
||||
//// StorageGradHessVector[j]*= -1.0;
|
||||
//// }
|
||||
//
|
||||
// Tensor<complex<double>,OHMMS_DIM> PG; PG=PrimLattice.G;
|
||||
// Tensor<complex<double>,OHMMS_DIM> TPG; TPG=transpose(PrimLattice.G);
|
||||
// Tensor<complex<double>,OHMMS_DIM> PG2;
|
||||
// // for (int n=0; n<OHMMS_DIM; n++) for (int m=0; m<OHMMS_DIM; m++) PG2(n,m)=TPG(n,m)*TPG(n,m);
|
||||
// PG2=dot(TPG,TPG);
|
||||
// Tensor<complex<double>,OHMMS_DIM> tmphs;
|
||||
// TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
//
|
||||
//
|
||||
// for (int j=0; j<NumValenceOrbs; j++) {
|
||||
// gradVec[j] = dot (PG, gradVec[j]);
|
||||
//// FIX FIX FIX: store transpose(PrimLattice.G)
|
||||
// tmphs = dot(TPG,StorageHessVector[j]);
|
||||
// hessVec[j] = dot(tmphs,PG);
|
||||
////Is this right?
|
||||
// tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
// gradhessVec[j]=dot(tmpghs,PG);
|
||||
//// gradhessVec[j]=dot(StorageGradHessVector[j],PrimLattice.G);
|
||||
// }
|
||||
//
|
||||
// // Add e^-ikr phase to B-spline orbitals
|
||||
// for (int j=0; j<NumValenceOrbs; j++) {
|
||||
// complex<double> u = valVec[j];
|
||||
// TinyVector<complex<double>,OHMMS_DIM> gradu = gradVec[j];
|
||||
// tmphs = hessVec[j];
|
||||
// tmpghs = gradhessVec[j];
|
||||
// PosType k = kPoints[j];
|
||||
// TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
// for (int n=0; n<OHMMS_DIM; n++) ck[n] = k[n];
|
||||
// double s,c;
|
||||
// double phase = -dot(r, k);
|
||||
// sincos (phase, &s, &c);
|
||||
// complex<double> e_mikr (c,s);
|
||||
//
|
||||
// valVec[j] = e_mikr*u;
|
||||
// gradVec[j] = e_mikr*(-eye*u*ck + gradu);
|
||||
// hessVec[j] = e_mikr*(tmphs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck));
|
||||
////Is this right?
|
||||
// for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
// for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
// for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
// hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
//
|
||||
// gradhessVec[j]=e_mikr*(tmpghs +eye*u*outerdot(ck,ck,ck) - symouterdot(gradu,ck,ck) -hvdot );
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// // Finally, copy into output vectors
|
||||
// int psiIndex = 0;
|
||||
// int N = StorageValueVector.size();
|
||||
// if (need2blend) {
|
||||
//
|
||||
// APP_ABORT("need2blend not implemented with Hessian evaluation.\n");
|
||||
// }
|
||||
// else { // No blending needed
|
||||
// for (int j=0; j<N; j++) {
|
||||
//// complex<double> psi_val;
|
||||
//// TinyVector<complex<double>, OHMMS_DIM> psi_grad;
|
||||
//// psi_val = StorageValueVector[j];
|
||||
//// psi_grad = StorageGradVector[j];
|
||||
//// tmphs = StorageHessVector[j];
|
||||
//// tmpghs = StorageGradHessVector[j];
|
||||
////
|
||||
// psi(i,psiIndex)=real(StorageValueVector[j]);
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// dpsi(i,psiIndex)(n) = real(StorageGradVector[j](n));
|
||||
// for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
// grad_grad_psi(i,psiIndex)(n) = real(StorageHessVector[j](n));
|
||||
//////Is this right?
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
// grad_grad_grad_logdet(i,psiIndex)[n](m) = real(StorageGradHessVector[j][n](m));
|
||||
// psiIndex++;
|
||||
//// // if (psiIndex >= dpsi.cols()) {
|
||||
//// // cerr << "Error: out of bounds writing in EinsplineSet::evalate.\n"
|
||||
//// // << "psiIndex = " << psiIndex << " dpsi.cols() = " << dpsi.cols() << endl;
|
||||
//// // }
|
||||
// if (MakeTwoCopies[j]) {
|
||||
// psi(i,psiIndex)=imag(StorageValueVector[j]);
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// dpsi(i,psiIndex)(n) = imag(StorageGradVector[j](n));
|
||||
// for (int n=0; n<OHMMS_DIM*OHMMS_DIM; n++)
|
||||
// grad_grad_psi(i,psiIndex)(n) = imag(StorageGradVector[j](n));
|
||||
//////Is this right?
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// for (int m=0; m<OHMMS_DIM*OHMMS_DIM; m++)
|
||||
// grad_grad_grad_logdet(i,psiIndex)[n](m) = imag(StorageGradHessVector[j][n](m));
|
||||
// psiIndex++;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// VGLMatTimer.stop();
|
||||
// }
|
||||
|
||||
template<> void
|
||||
EinsplineSetExtended<double>::evaluate_notranspose(const ParticleSet& P, int first, int last,
|
||||
ComplexValueMatrix_t& psi, ComplexGradMatrix_t& dpsi,
|
||||
ComplexHessMatrix_t& grad_grad_psi,
|
||||
ComplexGGGMatrix_t& grad_grad_grad_logdet)
|
||||
{
|
||||
APP_ABORT(" EinsplineSetExtended<StorageType>::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n");
|
||||
// VGLMatTimer.start();
|
||||
// for(int iat=first,i=0; iat<last; iat++,i++) {
|
||||
// PosType r (P.R[iat]);
|
||||
// PosType ru(PrimLattice.toUnit(P.R[iat]));
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// ru[n] -= std::floor (ru[n]);
|
||||
// EinsplineTimer.start();
|
||||
// EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
|
||||
// StorageGradVector, StorageHessVector,StorageGradHessVector);
|
||||
// EinsplineTimer.stop();
|
||||
// Tensor<double,OHMMS_DIM> TPG(transpose(PrimLattice.G));
|
||||
// Tensor<double,OHMMS_DIM> PG2;
|
||||
//// for (int n=0; n<OHMMS_DIM; n++) for (int m=0; m<OHMMS_DIM; m++) PG2(n,m)=TPG(n,m)*TPG(n,m);
|
||||
// PG2=dot(TPG,TPG);
|
||||
// TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
// complex<double> eye (0.0, 1.0);
|
||||
// for (int j=0; j<OrbitalSetSize; j++) {
|
||||
// complex<double> u;
|
||||
// TinyVector<complex<double>, OHMMS_DIM> gradu;
|
||||
// Tensor<complex<double>,OHMMS_DIM> hs,tmphs;
|
||||
//
|
||||
// u = StorageValueVector[j];
|
||||
// gradu = dot(PrimLattice.G, StorageGradVector[j]);
|
||||
// tmphs = dot(transpose(PrimLattice.G),StorageHessVector[j]);
|
||||
// hs = dot(tmphs,PrimLattice.G);
|
||||
// tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
// tmpghs = dot(tmpghs,PrimLattice.G);
|
||||
//
|
||||
// PosType k = kPoints[j];
|
||||
// TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
// for (int n=0; n<OHMMS_DIM; n++)
|
||||
// ck[n] = k[n];
|
||||
// double s,c;
|
||||
// double phase = -dot(r, k);
|
||||
// sincos (phase, &s, &c);
|
||||
// complex<double> e_mikr (c,s);
|
||||
// convert(e_mikr * u, psi(i,j));
|
||||
// //convert(e_mikr * u, psi(j,i));
|
||||
// convert(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
|
||||
// //convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
|
||||
// convert(e_mikr*(hs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck)),grad_grad_psi(i,j));
|
||||
// for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
// for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
// for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
// hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
//
|
||||
// convert(e_mikr*(tmpghs +eye*u*outerdot(ck,ck,ck)
|
||||
// -symouterdot(gradu,ck,ck)-hvdot ), grad_grad_grad_logdet(i,j));
|
||||
// }
|
||||
// }
|
||||
// VGLMatTimer.stop();
|
||||
}
|
||||
|
||||
template<typename StorageType> void
|
||||
EinsplineSetExtended<StorageType>::evaluate_notranspose(const ParticleSet& P, int first, int last,
|
||||
ComplexValueMatrix_t& psi, ComplexGradMatrix_t& dpsi,
|
||||
|
|
@ -1974,28 +2093,54 @@ namespace qmcplusplus {
|
|||
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
|
||||
StorageGradVector, StorageHessVector,StorageGradHessVector);
|
||||
EinsplineTimer.stop();
|
||||
Tensor<double,OHMMS_DIM> TPG(transpose(PrimLattice.G));
|
||||
Tensor<double,OHMMS_DIM> PG2;
|
||||
Tensor<complex<double>,OHMMS_DIM> PG; PG=PrimLattice.G;
|
||||
Tensor<complex<double>,OHMMS_DIM> TPG; TPG=transpose(PrimLattice.G);
|
||||
Tensor<complex<double>,OHMMS_DIM> PG2;
|
||||
// for (int n=0; n<OHMMS_DIM; n++) for (int m=0; m<OHMMS_DIM; m++) PG2(n,m)=TPG(n,m)*TPG(n,m);
|
||||
PG2=dot(TPG,TPG);
|
||||
Tensor<complex<double>,OHMMS_DIM> hs,tmphs;
|
||||
TinyVector<Tensor<complex<double>,OHMMS_DIM>,OHMMS_DIM> tmpghs,hvdot;
|
||||
for (int j=0; j<NumValenceOrbs; j++) {
|
||||
StorageGradVector[j] = dot (PG, StorageGradVector[j]);
|
||||
tmphs = dot(TPG,StorageHessVector[j]);
|
||||
StorageHessVector[j] = dot(tmphs,PG);
|
||||
tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
StorageGradHessVector[j]=dot(tmpghs,PG);
|
||||
}
|
||||
|
||||
complex<double> eye (0.0, 1.0);
|
||||
for (int j=0; j<OrbitalSetSize; j++) {
|
||||
complex<double> u;
|
||||
TinyVector<complex<double>, OHMMS_DIM> gradu;
|
||||
Tensor<complex<double>,OHMMS_DIM> hs,tmphs;
|
||||
complex<double> u = StorageValueVector[j];
|
||||
TinyVector<complex<double>,OHMMS_DIM> gradu = StorageGradVector[j];
|
||||
tmphs = StorageHessVector[j];
|
||||
tmpghs = StorageGradHessVector[j];
|
||||
|
||||
u = StorageValueVector[j];
|
||||
gradu = dot(PrimLattice.G, StorageGradVector[j]);
|
||||
tmphs = dot(transpose(PrimLattice.G),StorageHessVector[j]);
|
||||
hs = dot(tmphs,PrimLattice.G);
|
||||
tmpghs = dot(PG2,StorageGradHessVector[j]);
|
||||
tmpghs = dot(tmpghs,PrimLattice.G);
|
||||
PosType k = kPoints[j];
|
||||
TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
for (int n=0; n<OHMMS_DIM; n++) ck[n] = k[n];
|
||||
double s,c;
|
||||
double phase = -dot(r, k);
|
||||
sincos (phase, &s, &c);
|
||||
complex<double> e_mikr (c,s);
|
||||
|
||||
psi(i,j)=u;
|
||||
dpsi(i,j)=gradu;
|
||||
grad_grad_psi(i,j)=hs;
|
||||
grad_grad_grad_logdet(i,j)=tmpghs;
|
||||
StorageValueVector[j] = e_mikr*u;
|
||||
StorageGradVector[j] = e_mikr*(-eye*u*ck + gradu);
|
||||
StorageHessVector[j] = e_mikr*(tmphs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck));
|
||||
//Is this right?
|
||||
for(unsigned a0(0);a0<OHMMS_DIM;a0++)
|
||||
for(unsigned a1(0);a1<OHMMS_DIM;a1++)
|
||||
for(unsigned a2(0);a2<OHMMS_DIM;a2++)
|
||||
hvdot[a0](a1,a2) = eye*(tmphs(a0,a1)*ck[a2] + tmphs(a1,a2)*ck[a0] + tmphs(a0,a2)*ck[a1]);
|
||||
|
||||
|
||||
StorageGradHessVector[j]=e_mikr*(tmpghs +eye*u*outerdot(ck,ck,ck)
|
||||
-symouterdot(gradu,ck,ck) -hvdot );
|
||||
|
||||
|
||||
psi(i,j)=StorageValueVector[j];
|
||||
dpsi(i,j)=StorageGradVector[j];
|
||||
grad_grad_psi(i,j)=StorageHessVector[j];
|
||||
grad_grad_grad_logdet(i,j)=StorageGradHessVector[j];
|
||||
|
||||
// PosType k = kPoints[j];
|
||||
// TinyVector<complex<double>,OHMMS_DIM> ck;
|
||||
|
|
|
|||
|
|
@ -650,9 +650,9 @@ namespace qmcplusplus {
|
|||
temp=0.0;
|
||||
temp2=0.0;
|
||||
for(int j=0; j<NumPtcls; j++) {
|
||||
// for(int k=0; k<OHMMS_DIM; k++) temp2 += BFTrans->Bmat_full(i,FirstIndex+j)[k]*Fmat(j,j)[k];
|
||||
for(int k=0; k<OHMMS_DIM; k++) temp2 += BFTrans->Bmat_full(i,FirstIndex+j)[k]*Fmat(j,j)[k];
|
||||
temp += dot(BFTrans->Amat(i,FirstIndex+j),Fmat(j,j));
|
||||
temp2 += rcdot(BFTrans->Bmat_full(i,FirstIndex+j),Fmat(j,j));
|
||||
// temp2 += rcdot(BFTrans->Bmat_full(i,FirstIndex+j),Fmat(j,j));
|
||||
}
|
||||
myG(i) += temp;
|
||||
myL(i) += temp2;
|
||||
|
|
@ -811,7 +811,6 @@ namespace qmcplusplus {
|
|||
*/
|
||||
// Implementing the comments above
|
||||
Phi->evaluate(BFTrans->QP, FirstIndex, LastIndex, psiM, dpsiM, grad_grad_psiM, grad_grad_grad_psiM);
|
||||
|
||||
//std::copy(psiM.begin(),psiM.end(),psiMinv.begin());
|
||||
psiMinv=psiM;
|
||||
// invert backflow matrix
|
||||
|
|
@ -861,7 +860,7 @@ namespace qmcplusplus {
|
|||
PosType& cj = BFTrans->Cmat(pa,FirstIndex+j);
|
||||
for(int k=0; k<NumPtcls; k++) {
|
||||
f_a += (psiMinv(i,k)*dot(grad_grad_psiM(j,k),cj)
|
||||
- Fmat(k,j)*rcdot(BFTrans->Cmat(pa,FirstIndex+k),Fmat(i,k)));
|
||||
- Fmat(k,j)*dot(BFTrans->Cmat(pa,FirstIndex+k),Fmat(i,k)));
|
||||
}
|
||||
dFa(i,j)=f_a;
|
||||
}
|
||||
|
|
@ -875,9 +874,9 @@ namespace qmcplusplus {
|
|||
for(int j=0; j<NumPtcls; j++) {
|
||||
GradType B_j;
|
||||
for(int i=0; i<num; i++) B_j += BFTrans->Bmat_full(i,FirstIndex+j);
|
||||
dLa += (rcdot(Fmat(j,j),BFTrans->Ymat(pa,FirstIndex+j)) +
|
||||
dLa += (dot(Fmat(j,j),BFTrans->Ymat(pa,FirstIndex+j)) +
|
||||
dot(B_j,dFa(j,j)));
|
||||
dpsia += rcdot(Fmat(j,j),BFTrans->Cmat(pa,FirstIndex+j));
|
||||
dpsia += dot(Fmat(j,j),BFTrans->Cmat(pa,FirstIndex+j));
|
||||
}
|
||||
for(int j=0; j<NumPtcls; j++) {
|
||||
|
||||
|
|
@ -894,7 +893,7 @@ namespace qmcplusplus {
|
|||
q_j_prime += ( psiMinv(j,k)*(cj[0]*grad_grad_grad_psiM(j,k)[0]
|
||||
+ cj[1]*grad_grad_grad_psiM(j,k)[1]
|
||||
+ cj[2]*grad_grad_grad_psiM(j,k)[2])
|
||||
- rcdot(BFTrans->Cmat(pa,FirstIndex+k),Fmat(j,k))
|
||||
- dot(BFTrans->Cmat(pa,FirstIndex+k),Fmat(j,k))
|
||||
*Qmat(k,j) );
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -298,6 +298,7 @@ namespace qmcplusplus
|
|||
BackflowBuilder* bfbuilder = new BackflowBuilder(targetPtcl,ptclPool,targetPsi);
|
||||
bfbuilder->put(BFnode);
|
||||
BFTrans = bfbuilder->getBFTrans();
|
||||
delete bfbuilder;
|
||||
if(multiDet) {
|
||||
if(FastMSD) {
|
||||
multislaterdetfast_0->setBF(BFTrans);
|
||||
|
|
|
|||
|
|
@ -51,6 +51,7 @@ namespace qmcplusplus {
|
|||
void SPOSetBase::evaluate(const ParticleSet& P, int first, int last,
|
||||
ValueMatrix_t& logdet, GradMatrix_t& dlogdet, HessMatrix_t& grad_grad_logdet, GGGMatrix_t& grad_grad_grad_logdet)
|
||||
{
|
||||
logdet=0;
|
||||
evaluate_notranspose(P,first,last,t_logpsi,dlogdet,grad_grad_logdet,grad_grad_grad_logdet);
|
||||
transpose(t_logpsi.data(),logdet.data(),OrbitalSetSize);
|
||||
//evaluate_notranspose(P,first,last,logdet,dlogdet,d2logdet);
|
||||
|
|
|
|||
Loading…
Reference in New Issue