mirror of https://github.com/QMCPACK/qmcpack.git
Merge branch 'fix_opt_analyzer' of https://github.com/jtkrogel/qmcpack into fix_opt_analyzer
This commit is contained in:
Jaron Krogel
commit
021edb557c
|
|
@ -577,15 +577,11 @@ struct SplineC2CSoA: public SplineAdoptorBase<ST,3>
|
|||
template<typename VV, typename GV, typename GGV>
|
||||
void assign_vgh(const PointType& r, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
#if 0
|
||||
typedef std::complex<TT> ComplexT;
|
||||
CONSTEXPR ST zero(0);
|
||||
CONSTEXPR ST two(2);
|
||||
const ST g00=PrimLattice.G(0), g01=PrimLattice.G(1), g02=PrimLattice.G(2),
|
||||
g10=PrimLattice.G(3), g11=PrimLattice.G(4), g12=PrimLattice.G(5),
|
||||
g20=PrimLattice.G(6), g21=PrimLattice.G(7), g22=PrimLattice.G(8);
|
||||
const ST x=r[0], y=r[1], z=r[2];
|
||||
const ST symGG[6]={GGt[0],GGt[1]+GGt[3],GGt[2]+GGt[6],GGt[4],GGt[5]+GGt[7],GGt[8]};
|
||||
|
||||
const ST* restrict k0=myKcart.data(0);
|
||||
const ST* restrict k1=myKcart.data(1);
|
||||
|
|
@ -602,23 +598,8 @@ struct SplineC2CSoA: public SplineAdoptorBase<ST,3>
|
|||
const ST* restrict h22=myH.data(5);
|
||||
|
||||
const size_t N=kPoints.size();
|
||||
const size_t nsplines=myL.size();
|
||||
|
||||
ComplexT* restrict psi =psi.data(0)+first_spo; ASSUME_ALIGNED(psi);
|
||||
ComplexT* restrict vg_x=dpsi.data(0)+first_spo; ASSUME_ALIGNED(vg_x);
|
||||
ComplexT* restrict vg_y=dpsi.data(1)+first_spo; ASSUME_ALIGNED(vg_y);
|
||||
ComplexT* restrict vg_z=dpsi.data(2)+first_spo; ASSUME_ALIGNED(vg_z);
|
||||
ComplexT* restrict gg_xx=grad_grad_psi.data(0)+first_spo; ASSUME_ALIGNED(gg_xx);
|
||||
ComplexT* restrict gg_xy=grad_grad_psi.data(1)+first_spo; ASSUME_ALIGNED(gg_xy);
|
||||
ComplexT* restrict gg_xz=grad_grad_psi.data(2)+first_spo; ASSUME_ALIGNED(gg_xz);
|
||||
ComplexT* restrict gg_yx=grad_grad_psi.data(3)+first_spo; ASSUME_ALIGNED(gg_yx);
|
||||
ComplexT* restrict gg_yy=grad_grad_psi.data(4)+first_spo; ASSUME_ALIGNED(gg_yy);
|
||||
ComplexT* restrict gg_yz=grad_grad_psi.data(5)+first_spo; ASSUME_ALIGNED(gg_yz);
|
||||
ComplexT* restrict gg_zx=grad_grad_psi.data(6)+first_spo; ASSUME_ALIGNED(gg_zx);
|
||||
ComplexT* restrict gg_zy=grad_grad_psi.data(7)+first_spo; ASSUME_ALIGNED(gg_zy);
|
||||
ComplexT* restrict gg_zz=grad_grad_psi.data(8)+first_spo; ASSUME_ALIGNED(gg_zz);
|
||||
|
||||
#pragma simd
|
||||
#pragma omp simd
|
||||
for (size_t j=0; j<N; ++j)
|
||||
{
|
||||
int jr=j<<1;
|
||||
|
|
@ -651,42 +632,42 @@ struct SplineC2CSoA: public SplineAdoptorBase<ST,3>
|
|||
const ST gY_i=dY_i-val_r*kY;
|
||||
const ST gZ_i=dZ_i-val_r*kZ;
|
||||
|
||||
psi[j] =ComplexT(c*val_r-s*val_i,c*val_i+s*val_r);
|
||||
vg_x[j]=ComplexT(c*gX_r -s*gX_i, c*gX_i +s*gX_r);
|
||||
vg_y[j]=ComplexT(c*gY_r -s*gY_i, c*gY_i +s*gY_r);
|
||||
vg_z[j]=ComplexT(c*gZ_r -s*gZ_i, c*gZ_i +s*gZ_r);
|
||||
const size_t psiIndex=j+first_spo;
|
||||
psi[psiIndex] =ComplexT(c*val_r-s*val_i,c*val_i+s*val_r);
|
||||
dpsi[psiIndex][0]=ComplexT(c*gX_r -s*gX_i, c*gX_i +s*gX_r);
|
||||
dpsi[psiIndex][1]=ComplexT(c*gY_r -s*gY_i, c*gY_i +s*gY_r);
|
||||
dpsi[psiIndex][2]=ComplexT(c*gZ_r -s*gZ_i, c*gZ_i +s*gZ_r);
|
||||
|
||||
const ST kkV_r=mKK[j]*val_r; //kk*Vr
|
||||
const ST kkV_i=mKK[j]*val_i; //kk*Vi
|
||||
const ST h_xx_r=h00[jr]+kkV_r+kX*gX_r;
|
||||
const ST h_xy_r=h01[jr]+kkV_r+kX*gY_r;
|
||||
const ST h_xz_r=h02[jr]+kkV_r+kX*gZ_r;
|
||||
const ST h_yx_r=h01[jr]+kkV_r+kY*gX_r;
|
||||
const ST h_yy_r=h11[jr]+kkV_r+kY*gY_r;
|
||||
const ST h_yz_r=h12[jr]+kkV_r+kY*gZ_r;
|
||||
const ST h_zx_r=h02[jr]+kkV_r+kz*gX_r;
|
||||
const ST h_zy_r=h12[jr]+kkV_r+kz*gY_r;
|
||||
const ST h_zz_r=h22[jr]+kkV_r+kz*gZ_r;
|
||||
const ST h_xy_i=h01[ji]+kkV_i+kX*gY_i;
|
||||
const ST h_xz_i=h02[ji]+kkV_i+kX*gZ_i;
|
||||
const ST h_yx_i=h01[ji]+kkV_i+kY*gX_i;
|
||||
const ST h_yy_i=h11[ji]+kkV_i+kY*gY_i;
|
||||
const ST h_yz_i=h12[ji]+kkV_i+kY*gZ_i;
|
||||
const ST h_zx_i=h02[ji]+kkV_i+kz*gX_i;
|
||||
const ST h_zy_i=h12[ji]+kkV_i+kz*gY_i;
|
||||
const ST h_zz_i=h22[ji]+kkV_i+kz*gZ_i;
|
||||
const ST h_xx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g00,g01,g02)+kX*(gX_i+dX_i);
|
||||
const ST h_xy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g10,g11,g12)+kX*(gY_i+dY_i);
|
||||
const ST h_xz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g20,g21,g22)+kX*(gZ_i+dZ_i);
|
||||
const ST h_yx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g00,g01,g02)+kY*(gX_i+dX_i);
|
||||
const ST h_yy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g10,g11,g12)+kY*(gY_i+dY_i);
|
||||
const ST h_yz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g20,g21,g22)+kY*(gZ_i+dZ_i);
|
||||
const ST h_zx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g00,g01,g02)+kZ*(gX_i+dX_i);
|
||||
const ST h_zy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g10,g11,g12)+kZ*(gY_i+dY_i);
|
||||
const ST h_zz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g20,g21,g22)+kZ*(gZ_i+dZ_i);
|
||||
|
||||
gg_xx[j]=ComplexT(c*h_xx_r-s*h_xx_i, c*h_xx_i+s*h_xx_r);
|
||||
gg_xy[j]=ComplexT(c*h_xy_r-s*h_xy_i, c*h_xy_i+s*h_xy_r);
|
||||
gg_xz[j]=ComplexT(c*h_xz_r-s*h_xz_i, c*h_xz_i+s*h_xz_r);
|
||||
gg_yx[j]=ComplexT(c*h_yx_r-s*h_yx_i, c*h_yx_i+s*h_yx_r);
|
||||
gg_yy[j]=ComplexT(c*h_yy_r-s*h_yy_i, c*h_yy_i+s*h_yy_r);
|
||||
gg_yz[j]=ComplexT(c*h_yz_r-s*h_yz_i, c*h_yz_i+s*h_yz_r);
|
||||
gg_zx[j]=ComplexT(c*h_zx_r-s*h_zx_i, c*h_zx_i+s*h_zx_r);
|
||||
gg_zy[j]=ComplexT(c*h_zy_r-s*h_zy_i, c*h_zy_i+s*h_zy_r);
|
||||
gg_zz[j]=ComplexT(c*h_zz_r-s*h_zz_i, c*h_zz_i+s*h_zz_r);
|
||||
const ST h_xx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g00,g01,g02)-kX*(gX_r+dX_r);
|
||||
const ST h_xy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g10,g11,g12)-kX*(gY_r+dY_r);
|
||||
const ST h_xz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g20,g21,g22)-kX*(gZ_r+dZ_r);
|
||||
const ST h_yx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g00,g01,g02)-kY*(gX_r+dX_r);
|
||||
const ST h_yy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g10,g11,g12)-kY*(gY_r+dY_r);
|
||||
const ST h_yz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g20,g21,g22)-kY*(gZ_r+dZ_r);
|
||||
const ST h_zx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g00,g01,g02)-kZ*(gX_r+dX_r);
|
||||
const ST h_zy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g10,g11,g12)-kZ*(gY_r+dY_r);
|
||||
const ST h_zz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g20,g21,g22)-kZ*(gZ_r+dZ_r);
|
||||
|
||||
grad_grad_psi[psiIndex][0]=ComplexT(c*h_xx_r-s*h_xx_i, c*h_xx_i+s*h_xx_r);
|
||||
grad_grad_psi[psiIndex][1]=ComplexT(c*h_xy_r-s*h_xy_i, c*h_xy_i+s*h_xy_r);
|
||||
grad_grad_psi[psiIndex][2]=ComplexT(c*h_xz_r-s*h_xz_i, c*h_xz_i+s*h_xz_r);
|
||||
grad_grad_psi[psiIndex][3]=ComplexT(c*h_yx_r-s*h_yx_i, c*h_yx_i+s*h_yx_r);
|
||||
grad_grad_psi[psiIndex][4]=ComplexT(c*h_yy_r-s*h_yy_i, c*h_yy_i+s*h_yy_r);
|
||||
grad_grad_psi[psiIndex][5]=ComplexT(c*h_yz_r-s*h_yz_i, c*h_yz_i+s*h_yz_r);
|
||||
grad_grad_psi[psiIndex][6]=ComplexT(c*h_zx_r-s*h_zx_i, c*h_zx_i+s*h_zx_r);
|
||||
grad_grad_psi[psiIndex][7]=ComplexT(c*h_zy_r-s*h_zy_i, c*h_zy_i+s*h_zy_r);
|
||||
grad_grad_psi[psiIndex][8]=ComplexT(c*h_zz_r-s*h_zz_i, c*h_zz_i+s*h_zz_r);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
|
|
|
|||
|
|
@ -815,16 +815,10 @@ struct SplineC2RSoA: public SplineAdoptorBase<ST,3>
|
|||
template<typename VV, typename GV, typename GGV>
|
||||
void assign_vgh(const PointType& r, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
|
||||
#if 0
|
||||
typedef std::complex<TT> ComplexT;
|
||||
CONSTEXPR ST zero(0);
|
||||
CONSTEXPR ST two(2);
|
||||
const ST g00=PrimLattice.G(0), g01=PrimLattice.G(1), g02=PrimLattice.G(2),
|
||||
g10=PrimLattice.G(3), g11=PrimLattice.G(4), g12=PrimLattice.G(5),
|
||||
g20=PrimLattice.G(6), g21=PrimLattice.G(7), g22=PrimLattice.G(8);
|
||||
const ST x=r[0], y=r[1], z=r[2];
|
||||
const ST symGG[6]={GGt[0],GGt[1]+GGt[3],GGt[2]+GGt[6],GGt[4],GGt[5]+GGt[7],GGt[8]};
|
||||
|
||||
const ST* restrict k0=myKcart.data(0);
|
||||
const ST* restrict k1=myKcart.data(1);
|
||||
|
|
@ -841,24 +835,9 @@ struct SplineC2RSoA: public SplineAdoptorBase<ST,3>
|
|||
const ST* restrict h22=myH.data(5);
|
||||
|
||||
const size_t N=kPoints.size();
|
||||
const size_t nsplines=myL.size();
|
||||
|
||||
ComplexT* restrict psi =psi.data(0)+first_spo; ASSUME_ALIGNED(psi);
|
||||
ComplexT* restrict vg_x=dpsi.data(0)+first_spo; ASSUME_ALIGNED(vg_x);
|
||||
ComplexT* restrict vg_y=dpsi.data(1)+first_spo; ASSUME_ALIGNED(vg_y);
|
||||
ComplexT* restrict vg_z=dpsi.data(2)+first_spo; ASSUME_ALIGNED(vg_z);
|
||||
ComplexT* restrict gg_xx=grad_grad_psi.data(0)+first_spo; ASSUME_ALIGNED(gg_xx);
|
||||
ComplexT* restrict gg_xy=grad_grad_psi.data(1)+first_spo; ASSUME_ALIGNED(gg_xy);
|
||||
ComplexT* restrict gg_xz=grad_grad_psi.data(2)+first_spo; ASSUME_ALIGNED(gg_xz);
|
||||
ComplexT* restrict gg_yx=grad_grad_psi.data(3)+first_spo; ASSUME_ALIGNED(gg_yx);
|
||||
ComplexT* restrict gg_yy=grad_grad_psi.data(4)+first_spo; ASSUME_ALIGNED(gg_yy);
|
||||
ComplexT* restrict gg_yz=grad_grad_psi.data(5)+first_spo; ASSUME_ALIGNED(gg_yz);
|
||||
ComplexT* restrict gg_zx=grad_grad_psi.data(6)+first_spo; ASSUME_ALIGNED(gg_zx);
|
||||
ComplexT* restrict gg_zy=grad_grad_psi.data(7)+first_spo; ASSUME_ALIGNED(gg_zy);
|
||||
ComplexT* restrict gg_zz=grad_grad_psi.data(8)+first_spo; ASSUME_ALIGNED(gg_zz);
|
||||
|
||||
#pragma simd
|
||||
for (size_t j=0; j<N; ++j)
|
||||
#pragma omp simd
|
||||
for (size_t j=0; j<nComplexBands; j++)
|
||||
{
|
||||
int jr=j<<1;
|
||||
int ji=jr+1;
|
||||
|
|
@ -890,55 +869,129 @@ struct SplineC2RSoA: public SplineAdoptorBase<ST,3>
|
|||
const ST gY_i=dY_i-val_r*kY;
|
||||
const ST gZ_i=dZ_i-val_r*kZ;
|
||||
|
||||
psi[j] =ComplexT(c*val_r-s*val_i,c*val_i+s*val_r);
|
||||
vg_x[j]=ComplexT(c*gX_r -s*gX_i, c*gX_i +s*gX_r);
|
||||
vg_y[j]=ComplexT(c*gY_r -s*gY_i, c*gY_i +s*gY_r);
|
||||
vg_z[j]=ComplexT(c*gZ_r -s*gZ_i, c*gZ_i +s*gZ_r);
|
||||
const size_t psiIndex=first_spo+jr;
|
||||
|
||||
const ST kkV_r=mKK[j]*val_r; //kk*Vr
|
||||
const ST kkV_i=mKK[j]*val_i; //kk*Vi
|
||||
const ST h_xx_r=h00[jr]+kkV_r+kX*gX_r;
|
||||
const ST h_xy_r=h01[jr]+kkV_r+kX*gY_r;
|
||||
const ST h_xz_r=h02[jr]+kkV_r+kX*gZ_r;
|
||||
const ST h_yx_r=h01[jr]+kkV_r+kY*gX_r;
|
||||
const ST h_yy_r=h11[jr]+kkV_r+kY*gY_r;
|
||||
const ST h_yz_r=h12[jr]+kkV_r+kY*gZ_r;
|
||||
const ST h_zx_r=h02[jr]+kkV_r+kz*gX_r;
|
||||
const ST h_zy_r=h12[jr]+kkV_r+kz*gY_r;
|
||||
const ST h_zz_r=h22[jr]+kkV_r+kz*gZ_r;
|
||||
const ST h_xy_i=h01[ji]+kkV_i+kX*gY_i;
|
||||
const ST h_xz_i=h02[ji]+kkV_i+kX*gZ_i;
|
||||
const ST h_yx_i=h01[ji]+kkV_i+kY*gX_i;
|
||||
const ST h_yy_i=h11[ji]+kkV_i+kY*gY_i;
|
||||
const ST h_yz_i=h12[ji]+kkV_i+kY*gZ_i;
|
||||
const ST h_zx_i=h02[ji]+kkV_i+kz*gX_i;
|
||||
const ST h_zy_i=h12[ji]+kkV_i+kz*gY_i;
|
||||
const ST h_zz_i=h22[ji]+kkV_i+kz*gZ_i;
|
||||
psi[psiIndex] =c*val_r-s*val_i;
|
||||
dpsi[psiIndex][0] =c*gX_r-s*gX_i;
|
||||
dpsi[psiIndex][1] =c*gY_r-s*gY_i;
|
||||
dpsi[psiIndex][2] =c*gZ_r-s*gZ_i;
|
||||
|
||||
size_t psiIndex=jr; //only correct for j<nComplexbands, CORRECT ME
|
||||
gg_xx[psiIndex]=c*h_xx_r-s*h_xx_i;
|
||||
gg_xy[psiIndex]=c*h_xy_r-s*h_xy_i;
|
||||
gg_xz[psiIndex]=c*h_xz_r-s*h_xz_i;
|
||||
gg_yx[psiIndex]=c*h_yx_r-s*h_yx_i;
|
||||
gg_yy[psiIndex]=c*h_yy_r-s*h_yy_i;
|
||||
gg_yz[psiIndex]=c*h_yz_r-s*h_yz_i;
|
||||
gg_zx[psiIndex]=c*h_zx_r-s*h_zx_i;
|
||||
gg_zy[psiIndex]=c*h_zy_r-s*h_zy_i;
|
||||
gg_zz[psiIndex]=c*h_zz_r-s*h_zz_i;
|
||||
psi[psiIndex+1]=c*val_i+s*val_r;
|
||||
dpsi[psiIndex+1][0]=c*gX_i+s*gX_r;
|
||||
dpsi[psiIndex+1][1]=c*gY_i+s*gY_r;
|
||||
dpsi[psiIndex+1][2]=c*gZ_i+s*gZ_r;
|
||||
|
||||
if(j<nComplexBands) continue;
|
||||
const ST h_xx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g00,g01,g02)+kX*(gX_i+dX_i);
|
||||
const ST h_xy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g10,g11,g12)+kX*(gY_i+dY_i);
|
||||
const ST h_xz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g20,g21,g22)+kX*(gZ_i+dZ_i);
|
||||
const ST h_yx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g00,g01,g02)+kY*(gX_i+dX_i);
|
||||
const ST h_yy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g10,g11,g12)+kY*(gY_i+dY_i);
|
||||
const ST h_yz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g20,g21,g22)+kY*(gZ_i+dZ_i);
|
||||
const ST h_zx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g00,g01,g02)+kZ*(gX_i+dX_i);
|
||||
const ST h_zy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g10,g11,g12)+kZ*(gY_i+dY_i);
|
||||
const ST h_zz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g20,g21,g22)+kZ*(gZ_i+dZ_i);
|
||||
|
||||
gg_xx[psiIndex+1]=c*h_xx_i+s*h_xx_r;
|
||||
gg_xy[psiIndex+1]=c*h_xy_i+s*h_xy_r;
|
||||
gg_xz[psiIndex+1]=c*h_xz_i+s*h_xz_r;
|
||||
gg_yx[psiIndex+1]=c*h_yx_i+s*h_yx_r;
|
||||
gg_yy[psiIndex+1]=c*h_yy_i+s*h_yy_r;
|
||||
gg_yz[psiIndex+1]=c*h_yz_i+s*h_yz_r;
|
||||
gg_zx[psiIndex+1]=c*h_zx_i+s*h_zx_r;
|
||||
gg_zy[psiIndex+1]=c*h_zy_i+s*h_zy_r;
|
||||
gg_zz[psiIndex+1]=c*h_zz_i+s*h_zz_r;
|
||||
const ST h_xx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g00,g01,g02)-kX*(gX_r+dX_r);
|
||||
const ST h_xy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g10,g11,g12)-kX*(gY_r+dY_r);
|
||||
const ST h_xz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g20,g21,g22)-kX*(gZ_r+dZ_r);
|
||||
const ST h_yx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g00,g01,g02)-kY*(gX_r+dX_r);
|
||||
const ST h_yy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g10,g11,g12)-kY*(gY_r+dY_r);
|
||||
const ST h_yz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g20,g21,g22)-kY*(gZ_r+dZ_r);
|
||||
const ST h_zx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g00,g01,g02)-kZ*(gX_r+dX_r);
|
||||
const ST h_zy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g10,g11,g12)-kZ*(gY_r+dY_r);
|
||||
const ST h_zz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g20,g21,g22)-kZ*(gZ_r+dZ_r);
|
||||
|
||||
grad_grad_psi[psiIndex][0]=c*h_xx_r-s*h_xx_i;
|
||||
grad_grad_psi[psiIndex][1]=c*h_xy_r-s*h_xy_i;
|
||||
grad_grad_psi[psiIndex][2]=c*h_xz_r-s*h_xz_i;
|
||||
grad_grad_psi[psiIndex][3]=c*h_yx_r-s*h_yx_i;
|
||||
grad_grad_psi[psiIndex][4]=c*h_yy_r-s*h_yy_i;
|
||||
grad_grad_psi[psiIndex][5]=c*h_yz_r-s*h_yz_i;
|
||||
grad_grad_psi[psiIndex][6]=c*h_zx_r-s*h_zx_i;
|
||||
grad_grad_psi[psiIndex][7]=c*h_zy_r-s*h_zy_i;
|
||||
grad_grad_psi[psiIndex][8]=c*h_zz_r-s*h_zz_i;
|
||||
|
||||
grad_grad_psi[psiIndex+1][0]=c*h_xx_i+s*h_xx_r;
|
||||
grad_grad_psi[psiIndex+1][1]=c*h_xy_i+s*h_xy_r;
|
||||
grad_grad_psi[psiIndex+1][2]=c*h_xz_i+s*h_xz_r;
|
||||
grad_grad_psi[psiIndex+1][3]=c*h_yx_i+s*h_yx_r;
|
||||
grad_grad_psi[psiIndex+1][4]=c*h_yy_i+s*h_yy_r;
|
||||
grad_grad_psi[psiIndex+1][5]=c*h_yz_i+s*h_yz_r;
|
||||
grad_grad_psi[psiIndex+1][6]=c*h_zx_i+s*h_zx_r;
|
||||
grad_grad_psi[psiIndex+1][7]=c*h_zy_i+s*h_zy_r;
|
||||
grad_grad_psi[psiIndex+1][8]=c*h_zz_i+s*h_zz_r;
|
||||
}
|
||||
|
||||
#pragma omp simd
|
||||
for (size_t j=nComplexBands; j<N; j++)
|
||||
{
|
||||
int jr=j<<1;
|
||||
int ji=jr+1;
|
||||
|
||||
const ST kX=k0[j];
|
||||
const ST kY=k1[j];
|
||||
const ST kZ=k2[j];
|
||||
const ST val_r=myV[jr];
|
||||
const ST val_i=myV[ji];
|
||||
|
||||
//phase
|
||||
ST s, c;
|
||||
sincos(-(x*kX+y*kY+z*kZ),&s,&c);
|
||||
|
||||
//dot(PrimLattice.G,myG[j])
|
||||
const ST dX_r = g00*g0[jr]+g01*g1[jr]+g02*g2[jr];
|
||||
const ST dY_r = g10*g0[jr]+g11*g1[jr]+g12*g2[jr];
|
||||
const ST dZ_r = g20*g0[jr]+g21*g1[jr]+g22*g2[jr];
|
||||
|
||||
const ST dX_i = g00*g0[ji]+g01*g1[ji]+g02*g2[ji];
|
||||
const ST dY_i = g10*g0[ji]+g11*g1[ji]+g12*g2[ji];
|
||||
const ST dZ_i = g20*g0[ji]+g21*g1[ji]+g22*g2[ji];
|
||||
|
||||
// \f$\nabla \psi_r + {\bf k}\psi_i\f$
|
||||
const ST gX_r=dX_r+val_i*kX;
|
||||
const ST gY_r=dY_r+val_i*kY;
|
||||
const ST gZ_r=dZ_r+val_i*kZ;
|
||||
const ST gX_i=dX_i-val_r*kX;
|
||||
const ST gY_i=dY_i-val_r*kY;
|
||||
const ST gZ_i=dZ_i-val_r*kZ;
|
||||
|
||||
const size_t psiIndex=first_spo+nComplexBands+j;
|
||||
|
||||
psi[psiIndex] =c*val_r-s*val_i;
|
||||
dpsi[psiIndex][0] =c*gX_r-s*gX_i;
|
||||
dpsi[psiIndex][1] =c*gY_r-s*gY_i;
|
||||
dpsi[psiIndex][2] =c*gZ_r-s*gZ_i;
|
||||
|
||||
const ST h_xx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g00,g01,g02)+kX*(gX_i+dX_i);
|
||||
const ST h_xy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g10,g11,g12)+kX*(gY_i+dY_i);
|
||||
const ST h_xz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g00,g01,g02,g20,g21,g22)+kX*(gZ_i+dZ_i);
|
||||
const ST h_yx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g00,g01,g02)+kY*(gX_i+dX_i);
|
||||
const ST h_yy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g10,g11,g12)+kY*(gY_i+dY_i);
|
||||
const ST h_yz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g10,g11,g12,g20,g21,g22)+kY*(gZ_i+dZ_i);
|
||||
const ST h_zx_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g00,g01,g02)+kZ*(gX_i+dX_i);
|
||||
const ST h_zy_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g10,g11,g12)+kZ*(gY_i+dY_i);
|
||||
const ST h_zz_r=v_m_v(h00[jr],h01[jr],h02[jr],h11[jr],h12[jr],h22[jr],g20,g21,g22,g20,g21,g22)+kZ*(gZ_i+dZ_i);
|
||||
|
||||
const ST h_xx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g00,g01,g02)-kX*(gX_r+dX_r);
|
||||
const ST h_xy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g10,g11,g12)-kX*(gY_r+dY_r);
|
||||
const ST h_xz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g00,g01,g02,g20,g21,g22)-kX*(gZ_r+dZ_r);
|
||||
const ST h_yx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g00,g01,g02)-kY*(gX_r+dX_r);
|
||||
const ST h_yy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g10,g11,g12)-kY*(gY_r+dY_r);
|
||||
const ST h_yz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g10,g11,g12,g20,g21,g22)-kY*(gZ_r+dZ_r);
|
||||
const ST h_zx_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g00,g01,g02)-kZ*(gX_r+dX_r);
|
||||
const ST h_zy_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g10,g11,g12)-kZ*(gY_r+dY_r);
|
||||
const ST h_zz_i=v_m_v(h00[ji],h01[ji],h02[ji],h11[ji],h12[ji],h22[ji],g20,g21,g22,g20,g21,g22)-kZ*(gZ_r+dZ_r);
|
||||
|
||||
grad_grad_psi[psiIndex][0]=c*h_xx_r-s*h_xx_i;
|
||||
grad_grad_psi[psiIndex][1]=c*h_xy_r-s*h_xy_i;
|
||||
grad_grad_psi[psiIndex][2]=c*h_xz_r-s*h_xz_i;
|
||||
grad_grad_psi[psiIndex][3]=c*h_yx_r-s*h_yx_i;
|
||||
grad_grad_psi[psiIndex][4]=c*h_yy_r-s*h_yy_i;
|
||||
grad_grad_psi[psiIndex][5]=c*h_yz_r-s*h_yz_i;
|
||||
grad_grad_psi[psiIndex][6]=c*h_zx_r-s*h_zx_i;
|
||||
grad_grad_psi[psiIndex][7]=c*h_zy_r-s*h_zy_i;
|
||||
grad_grad_psi[psiIndex][8]=c*h_zz_r-s*h_zz_i;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
|
|
|
|||
|
|
@ -433,12 +433,9 @@ struct SplineR2RSoA: public SplineAdoptorBase<ST,3>
|
|||
void assign_vgh(int bc_sign, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
const ST cone = (bc_sign &1)? -1:1;
|
||||
const ST g00=cone*PrimLattice.G(0), g01=cone*PrimLattice.G(1), g02=cone*PrimLattice.G(2),
|
||||
g10=cone*PrimLattice.G(3), g11=cone*PrimLattice.G(4), g12=cone*PrimLattice.G(5),
|
||||
g20=cone*PrimLattice.G(6), g21=cone*PrimLattice.G(7), g22=cone*PrimLattice.G(8);
|
||||
const ST gg00=cone*GGt(0), gg01=cone*GGt(1), gg02=cone*GGt(2),
|
||||
gg10=cone*GGt(3), gg11=cone*GGt(4), gg12=cone*GGt(5),
|
||||
gg20=cone*GGt(6), gg21=cone*GGt(7), gg22=cone*GGt(8);
|
||||
const ST g00=PrimLattice.G(0), g01=PrimLattice.G(1), g02=PrimLattice.G(2),
|
||||
g10=PrimLattice.G(3), g11=PrimLattice.G(4), g12=PrimLattice.G(5),
|
||||
g20=PrimLattice.G(6), g21=PrimLattice.G(7), g22=PrimLattice.G(8);
|
||||
|
||||
const ST* restrict g0=myG.data(0);
|
||||
const ST* restrict g1=myG.data(1);
|
||||
|
|
@ -450,26 +447,41 @@ struct SplineR2RSoA: public SplineAdoptorBase<ST,3>
|
|||
const ST* restrict h12=myH.data(4);
|
||||
const ST* restrict h22=myH.data(5);
|
||||
|
||||
const size_t nsplines=myL.size();
|
||||
const size_t N=last_spo-first_spo;
|
||||
#pragma omp simd
|
||||
for(size_t j=0; j<N; ++j)
|
||||
|
||||
#pragma omp simd
|
||||
for (size_t j=0; j<N; ++j)
|
||||
{
|
||||
size_t psiIndex=first_spo+j;
|
||||
psi[psiIndex]=cone*myV[j];
|
||||
dpsi[psiIndex][0]=(g00*g0[j]+g01*g1[j]+g02*g2[j]);
|
||||
dpsi[psiIndex][1]=(g10*g0[j]+g11*g1[j]+g12*g2[j]);
|
||||
dpsi[psiIndex][2]=(g20*g0[j]+g21*g1[j]+g22*g2[j]);
|
||||
//grad_grad_psi[psiIndex]=dot(myH,GGt)
|
||||
grad_grad_psi[psiIndex][0]=h00[j]*gg00+h01[j]*gg10+h02[j]*gg20;
|
||||
grad_grad_psi[psiIndex][1]=h00[j]*gg01+h01[j]*gg11+h02[j]*gg21;
|
||||
grad_grad_psi[psiIndex][2]=h00[j]*gg02+h01[j]*gg12+h02[j]*gg22;
|
||||
grad_grad_psi[psiIndex][3]=h01[j]*gg00+h11[j]*gg10+h12[j]*gg20;
|
||||
grad_grad_psi[psiIndex][4]=h01[j]*gg01+h11[j]*gg11+h12[j]*gg21;
|
||||
grad_grad_psi[psiIndex][5]=h01[j]*gg02+h11[j]*gg12+h12[j]*gg22;
|
||||
grad_grad_psi[psiIndex][6]=h02[j]*gg00+h12[j]*gg10+h22[j]*gg20;
|
||||
grad_grad_psi[psiIndex][7]=h02[j]*gg01+h12[j]*gg11+h22[j]*gg21;
|
||||
grad_grad_psi[psiIndex][8]=h02[j]*gg02+h12[j]*gg12+h22[j]*gg22;
|
||||
//dot(PrimLattice.G,myG[j])
|
||||
const ST dX_r = g00*g0[j]+g01*g1[j]+g02*g2[j];
|
||||
const ST dY_r = g10*g0[j]+g11*g1[j]+g12*g2[j];
|
||||
const ST dZ_r = g20*g0[j]+g21*g1[j]+g22*g2[j];
|
||||
|
||||
const size_t psiIndex=j+first_spo;
|
||||
psi[psiIndex] =cone*myV[j];
|
||||
dpsi[psiIndex][0]=cone*dX_r;
|
||||
dpsi[psiIndex][1]=cone*dY_r;
|
||||
dpsi[psiIndex][2]=cone*dZ_r;
|
||||
|
||||
const ST h_xx_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g00,g01,g02,g00,g01,g02);
|
||||
const ST h_xy_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g00,g01,g02,g10,g11,g12);
|
||||
const ST h_xz_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g00,g01,g02,g20,g21,g22);
|
||||
const ST h_yx_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g10,g11,g12,g00,g01,g02);
|
||||
const ST h_yy_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g10,g11,g12,g10,g11,g12);
|
||||
const ST h_yz_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g10,g11,g12,g20,g21,g22);
|
||||
const ST h_zx_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g20,g21,g22,g00,g01,g02);
|
||||
const ST h_zy_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g20,g21,g22,g10,g11,g12);
|
||||
const ST h_zz_r=v_m_v(h00[j],h01[j],h02[j],h11[j],h12[j],h22[j],g20,g21,g22,g20,g21,g22);
|
||||
|
||||
grad_grad_psi[psiIndex][0]=cone*h_xx_r;
|
||||
grad_grad_psi[psiIndex][1]=cone*h_xy_r;
|
||||
grad_grad_psi[psiIndex][2]=cone*h_xz_r;
|
||||
grad_grad_psi[psiIndex][3]=cone*h_yx_r;
|
||||
grad_grad_psi[psiIndex][4]=cone*h_yy_r;
|
||||
grad_grad_psi[psiIndex][5]=cone*h_yz_r;
|
||||
grad_grad_psi[psiIndex][6]=cone*h_zx_r;
|
||||
grad_grad_psi[psiIndex][7]=cone*h_zy_r;
|
||||
grad_grad_psi[psiIndex][8]=cone*h_zz_r;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -9,19 +9,15 @@
|
|||
// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp.
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "QMCWaveFunctions/BsplineFactory/macro.h"
|
||||
#include "Numerics/e2iphi.h"
|
||||
#include "simd/vmath.hpp"
|
||||
#include "qmc_common.h"
|
||||
#include <Utilities/ProgressReportEngine.h>
|
||||
#include "QMCWaveFunctions/EinsplineSetBuilder.h"
|
||||
#include "QMCWaveFunctions/EinsplineAdoptor.h"
|
||||
#include "QMCWaveFunctions/SplineC2XAdoptor.h"
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineC2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineC2CAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/HybridCplxAdoptor.h"
|
||||
#endif
|
||||
#include <fftw3.h>
|
||||
#include <QMCWaveFunctions/einspline_helper.hpp>
|
||||
#include "QMCWaveFunctions/BsplineReaderBase.h"
|
||||
|
|
@ -37,25 +33,15 @@ namespace qmcplusplus
|
|||
BsplineReaderBase* aReader=nullptr;
|
||||
|
||||
#if defined(QMC_COMPLEX)
|
||||
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridCplxSoA<SplineC2CSoA<double,RealType> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineC2CSoA<double,RealType> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineC2CPackedAdoptor<double,RealType,3> >(e);
|
||||
#endif
|
||||
#else //QMC_COMPLEX
|
||||
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridCplxSoA<SplineC2RSoA<double,RealType> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineC2RSoA<double,RealType> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineC2RPackedAdoptor<double,RealType,3> >(e);
|
||||
#endif
|
||||
#endif
|
||||
|
||||
return aReader;
|
||||
|
|
|
|||
|
|
@ -8,19 +8,16 @@
|
|||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp.
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
#include "QMCWaveFunctions/BsplineFactory/macro.h"
|
||||
|
||||
#include "Numerics/e2iphi.h"
|
||||
#include "simd/vmath.hpp"
|
||||
#include "qmc_common.h"
|
||||
#include <Utilities/ProgressReportEngine.h>
|
||||
#include "QMCWaveFunctions/EinsplineSetBuilder.h"
|
||||
#include "QMCWaveFunctions/EinsplineAdoptor.h"
|
||||
#include "QMCWaveFunctions/SplineC2XAdoptor.h"
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineC2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineC2CAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/HybridCplxAdoptor.h"
|
||||
#endif
|
||||
#include <fftw3.h>
|
||||
#include <QMCWaveFunctions/einspline_helper.hpp>
|
||||
#include "QMCWaveFunctions/BsplineReaderBase.h"
|
||||
|
|
@ -36,24 +33,15 @@ namespace qmcplusplus
|
|||
BsplineReaderBase* aReader=nullptr;
|
||||
|
||||
#if defined(QMC_COMPLEX)
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridCplxSoA<SplineC2CSoA<float,RealType> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineC2CSoA<float,RealType> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineC2CPackedAdoptor<float,RealType,3> >(e);
|
||||
#endif
|
||||
#else //QMC_COMPLEX
|
||||
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridCplxSoA<SplineC2RSoA<float,RealType> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineC2RSoA<float,RealType> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineC2RPackedAdoptor<float,RealType,3> >(e);
|
||||
#endif
|
||||
#endif
|
||||
|
||||
return aReader;
|
||||
|
|
|
|||
|
|
@ -8,12 +8,11 @@
|
|||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp.
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
#include "QMCWaveFunctions/BsplineFactory/macro.h"
|
||||
|
||||
#include "qmc_common.h"
|
||||
#include <Utilities/ProgressReportEngine.h>
|
||||
#include "QMCWaveFunctions/EinsplineSetBuilder.h"
|
||||
#include "QMCWaveFunctions/EinsplineAdoptor.h"
|
||||
#include "QMCWaveFunctions/SplineR2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineR2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/HybridRealAdoptor.h"
|
||||
#include <fftw3.h>
|
||||
|
|
@ -28,14 +27,10 @@ namespace qmcplusplus
|
|||
BsplineReaderBase* createBsplineRealDouble(EinsplineSetBuilder* e, bool hybrid_rep)
|
||||
{
|
||||
BsplineReaderBase* aReader=nullptr;
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridRealSoA<SplineR2RSoA<double,OHMMS_PRECISION> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineR2RSoA<double,OHMMS_PRECISION> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineR2RAdoptor<double,OHMMS_PRECISION,3> >(e);
|
||||
#endif
|
||||
return aReader;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -8,12 +8,11 @@
|
|||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp.
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "qmc_common.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/macro.h"
|
||||
#include <Utilities/ProgressReportEngine.h>
|
||||
#include "QMCWaveFunctions/EinsplineSetBuilder.h"
|
||||
#include "QMCWaveFunctions/EinsplineAdoptor.h"
|
||||
#include "QMCWaveFunctions/SplineR2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/SplineR2RAdoptor.h"
|
||||
#include "QMCWaveFunctions/BsplineFactory/HybridRealAdoptor.h"
|
||||
#include <fftw3.h>
|
||||
|
|
@ -29,14 +28,10 @@ namespace qmcplusplus
|
|||
{
|
||||
BsplineReaderBase* aReader=nullptr;
|
||||
|
||||
#if defined(QMC_ENABLE_SOA_DET)
|
||||
if(hybrid_rep)
|
||||
aReader= new SplineHybridAdoptorReader<HybridRealSoA<SplineR2RSoA<float,OHMMS_PRECISION> > >(e);
|
||||
else
|
||||
aReader= new SplineAdoptorReader<SplineR2RSoA<float,OHMMS_PRECISION> >(e);
|
||||
#else
|
||||
aReader= new SplineAdoptorReader<SplineR2RAdoptor<float,OHMMS_PRECISION,3> >(e);
|
||||
#endif
|
||||
return aReader;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,3 +0,0 @@
|
|||
#ifndef QMC_ENABLE_SOA_DET
|
||||
#define QMC_ENABLE_SOA_DET 1
|
||||
#endif
|
||||
|
|
@ -297,8 +297,7 @@ EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur)
|
|||
kid = kid->next;
|
||||
}
|
||||
|
||||
if (has_backflow && UseRealOrbitals) APP_ABORT("backflow optimization is broken with UseRealOrbitals");
|
||||
if (has_backflow && use_einspline_set_extended!="yes") APP_ABORT("backflow optimization does not yet function with EinsplinAdoptor, please add use_old_spline=\"yes\" to <determinantset> or <sposet_builder>.");
|
||||
if (has_backflow && use_einspline_set_extended=="yes" && UseRealOrbitals) APP_ABORT("backflow optimization is broken with UseRealOrbitals");
|
||||
|
||||
//////////////////////////////////
|
||||
// Create the OrbitalSet object
|
||||
|
|
|
|||
|
|
@ -18,8 +18,8 @@
|
|||
*
|
||||
* The most general reader class for SplineAdoptor using the full single grid for the supercell
|
||||
* - SplineR2RAdoptor
|
||||
* - SplineC2CPackedAdoptor
|
||||
* - SplineC2RPackedAdoptor
|
||||
* - SplineC2CAdoptor
|
||||
* - SplineC2RAdoptor
|
||||
* Each band is initialized with UBspline_3d_d and both real and imaginary parts are passed to the adoptor object
|
||||
* which will convert the data type to their internal precision.
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -1,785 +0,0 @@
|
|||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
// This file is distributed under the University of Illinois/NCSA Open Source License.
|
||||
// See LICENSE file in top directory for details.
|
||||
//
|
||||
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
|
||||
//
|
||||
// File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Ye Luo, yeluo@anl.gov, Argonne National Laboratory
|
||||
// Anouar Benali, benali@anl.gov, Argonne National Laboratory
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
/** @file SplineC2XAdoptor.h
|
||||
*
|
||||
* Adoptor classes to handle complex-to-(real,complex) with arbitrary precision
|
||||
*/
|
||||
#ifndef QMCPLUSPLUS_EINSPLINE_C2X_ADOPTOR_H
|
||||
#define QMCPLUSPLUS_EINSPLINE_C2X_ADOPTOR_H
|
||||
|
||||
#include <spline/einspline_engine.hpp>
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
//template<typename T, unsigned D>
|
||||
// inline Tensor<T,D> dot(const Tensor<T,D>& a, const Tensor<T,D>& b)
|
||||
// {
|
||||
// Tensor<T,D> res;
|
||||
// for(int i=0; i<D; ++i)
|
||||
// for(int j=0; j<D; ++j)
|
||||
// {
|
||||
// T s=0;
|
||||
// for(int k=0; k<D; ++k)
|
||||
// s+=a(i,k)*b(k,j);
|
||||
// res(i,j)=s;
|
||||
// }
|
||||
// return res;
|
||||
// }
|
||||
|
||||
|
||||
/** adoptor class to match std::complex<ST> spline stored in a packed array with std::complex<TT> SPOs
|
||||
* @tparam ST precision of spline
|
||||
* @tparam TT precision of SPOs
|
||||
* @tparam D dimension
|
||||
*/
|
||||
template<typename ST, typename TT, unsigned D>
|
||||
struct SplineC2CPackedAdoptor: public SplineAdoptorBase<ST,D>
|
||||
{
|
||||
typedef typename einspline_traits<ST,D>::SplineType SplineType;
|
||||
typedef typename einspline_traits<ST,D>::BCType BCType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::PointType PointType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::SingleSplineType SingleSplineType;
|
||||
|
||||
using SplineAdoptorBase<ST,D>::first_spo;
|
||||
using SplineAdoptorBase<ST,D>::last_spo;
|
||||
using SplineAdoptorBase<ST,D>::GGt;
|
||||
using SplineAdoptorBase<ST,D>::PrimLattice;
|
||||
using SplineAdoptorBase<ST,D>::kPoints;
|
||||
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myV;
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myL;
|
||||
typename OrbitalSetTraits<ST>::GradVector_t myG;
|
||||
typename OrbitalSetTraits<ST>::HessVector_t myH;
|
||||
typename OrbitalSetTraits<ST>::GradHessVector_t myGH;
|
||||
|
||||
///Actual spline table, multi_bspline_3d_(d,s)
|
||||
SplineType *MultiSpline;
|
||||
///number of points of the original grid
|
||||
int BaseN[3];
|
||||
///offset of the original grid, always 0
|
||||
int BaseOffset[3];
|
||||
|
||||
//vector<ST> phase;
|
||||
|
||||
SplineC2CPackedAdoptor():MultiSpline(0)
|
||||
{
|
||||
this->is_complex=true;
|
||||
this->AdoptorName="SplineC2CPackedAdoptor";
|
||||
this->KeyWord="C2CPacked";
|
||||
}
|
||||
|
||||
inline void resizeStorage(int n, int nvals)
|
||||
{
|
||||
SplineAdoptorBase<ST,D>::init_base(n);
|
||||
myV.resize(2*n);
|
||||
myL.resize(2*n);
|
||||
myG.resize(2*n);
|
||||
myH.resize(2*n);
|
||||
}
|
||||
|
||||
template<typename GT, typename BCT>
|
||||
void create_spline(GT& xyz_g, BCT& xyz_bc)
|
||||
{
|
||||
MultiSpline=einspline::create(MultiSpline,xyz_g,xyz_bc,myV.size());
|
||||
for(int i=0; i<D; ++i)
|
||||
{
|
||||
BaseOffset[i]=0;
|
||||
BaseN[i]=xyz_g[i].num+3;
|
||||
}
|
||||
qmc_common.memory_allocated += MultiSpline->coefs_size*sizeof(ST);
|
||||
}
|
||||
|
||||
inline void set_spline(ST* restrict psi_r, ST* restrict psi_i, int twist, int ispline, int level)
|
||||
{
|
||||
einspline::set(MultiSpline, 2*ispline, psi_r);
|
||||
einspline::set(MultiSpline, 2*ispline+1, psi_i);
|
||||
}
|
||||
|
||||
inline void set_spline(SingleSplineType* spline_r, SingleSplineType* spline_i, int twist, int ispline, int level)
|
||||
{
|
||||
einspline::set(MultiSpline, 2*ispline, spline_r, BaseOffset, BaseN);
|
||||
einspline::set(MultiSpline, 2*ispline+1,spline_i, BaseOffset, BaseN);
|
||||
}
|
||||
|
||||
inline void set_spline_domain(SingleSplineType* spline_r, SingleSplineType* spline_i,
|
||||
int twist, int ispline, const int* offset_l, const int* mesh_l)
|
||||
{
|
||||
einspline::set(MultiSpline, 2*ispline, spline_r, offset_l, mesh_l);
|
||||
einspline::set(MultiSpline, 2*ispline+1,spline_i, offset_l, mesh_l);
|
||||
}
|
||||
|
||||
bool read_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.read(bigtable,o.str().c_str()); //"spline_0");
|
||||
}
|
||||
|
||||
bool write_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.write(bigtable,o.str().c_str()); //"spline_0");
|
||||
}
|
||||
|
||||
/** assign myV to psi
|
||||
*
|
||||
* Taken out for the derived classes
|
||||
*/
|
||||
template<typename VV>
|
||||
inline void assign_v(const PointType& r, int bc_sign, VV& psi)
|
||||
{
|
||||
register ST s,c;
|
||||
//TT* restrict t_ptr=reinterpret_cast<TT*>(psi.data())+(first_spo<<1);
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
{
|
||||
int jr=j<<1;
|
||||
sincos(-dot(r,kPoints[j]),&s,&c);
|
||||
psi[psiIndex]= std::complex<TT>(c*myV[jr]-s*myV[jr+1],s*myV[jr]+c*myV[jr+1]);
|
||||
//t_ptr[jr ]=c*myV[jr]-s*myV[jr+1];
|
||||
//t_ptr[jr+1]=s*myV[jr]+c*myV[jr+1];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV>
|
||||
inline void evaluate_v(const ParticleSet& P, const int iat, VV& psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate(MultiSpline,ru,myV);
|
||||
assign_v(r,0,psi);
|
||||
////computePhases(r);
|
||||
////simd::multiadd(psi.size(),eikr.data(),psi.data());
|
||||
}
|
||||
|
||||
/** assign internal data to psi,dpsi,d2psi
|
||||
*/
|
||||
template<typename VV, typename GV>
|
||||
inline void assign_vgl(const PointType& r, int bc_sign, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
const int N=kPoints.size();
|
||||
for (int j=0; j<2*N; j++)
|
||||
myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
for (int j=0; j<2*N; j++)
|
||||
myL[j] = trace(myH[j],GGt);
|
||||
const ST two=2.0;
|
||||
ST s,c;
|
||||
PointType g_r, g_i;
|
||||
//can easily make three independent loops
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
{
|
||||
int jr=j<<1;
|
||||
int ji=jr+1;
|
||||
g_r=myG[jr]+myV[ji]*kPoints[j]; // \f$\nabla \psi_r + {\bf k}\psi_i\f$
|
||||
g_i=myG[ji]-myV[jr]*kPoints[j]; // \f$\nabla \psi_i - {\bf k}\psi_r\f$
|
||||
ST kk=-dot(kPoints[j],kPoints[j]);
|
||||
myL[jr]+=kk*myV[jr]+two*dot(kPoints[j],myG[ji]);
|
||||
myL[ji]+=kk*myV[ji]-two*dot(kPoints[j],myG[jr]);
|
||||
sincos(-dot(r,kPoints[j]),&s,&c); //e-ikr (beware of -1)
|
||||
psi[psiIndex]= std::complex<TT>(c*myV[jr]-s*myV[ji],c*myV[ji]+s*myV[jr]);
|
||||
d2psi[psiIndex]= std::complex<TT>(c*myL[jr]-s*myL[ji],c*myL[ji]+s*myL[jr]);
|
||||
for(int idim=0; idim<D; ++idim)
|
||||
dpsi[psiIndex][idim]= std::complex<TT>(c*g_r[idim]-s*g_i[idim], c*g_i[idim]+s*g_r[idim]);
|
||||
}
|
||||
//complex<ST> e_mikr(c,s);
|
||||
//convert(e_mikr * myV[j], psi[j]);
|
||||
//convert(e_mikr*(-myV[j]*ck + myG[j]), dpsi[j]);
|
||||
//convert(e_mikr*(-myV[j]*kk - two*dot(ck,myG[j]) + myL[j]), d2psi[j]);
|
||||
}
|
||||
|
||||
template<typename VV, typename GV>
|
||||
inline void evaluate_vgl(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgl(r,0,psi,dpsi,d2psi);
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void assign_vgh(const PointType& r, int bc_sign, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
//convert to Cartesian G and Hess
|
||||
const int N=kPoints.size();
|
||||
for (int j=0; j<2*N; j++)
|
||||
myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
for (int j=0; j<2*N; j++)
|
||||
myH[j] = dot(myH[j],GGt);
|
||||
ST s,c;
|
||||
PointType g_r, g_i;
|
||||
Tensor<ST,D> kk,h_r,h_i;
|
||||
//can easily make three independent loops
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
{
|
||||
int jr=j<<1;
|
||||
int ji=jr+1;
|
||||
g_r=myG[jr]+myV[ji]*kPoints[j]; // \f$\nabla \psi_r + {\bf k}\psi_i\f$
|
||||
g_i=myG[ji]-myV[jr]*kPoints[j]; // \f$\nabla \psi_i - {\bf k}\psi_r\f$
|
||||
sincos(-dot(r,kPoints[j]),&s,&c); //e-ikr (beware of -1)
|
||||
psi[psiIndex]= std::complex<TT>(c*myV[jr]-s*myV[ji],c*myV[ji]+s*myV[jr]);
|
||||
for(int idim=0; idim<D; ++idim)
|
||||
dpsi[psiIndex][idim]= std::complex<TT>(c*g_r[idim]-s*g_i[idim], c*g_i[idim]+s*g_r[idim]);
|
||||
kk=outerProduct(kPoints[j],kPoints[j]); // \f$kk=k^k \f$
|
||||
h_r=myH[jr]-myV[jr]*kk+outerProductSymm(kPoints[j],myG[ji]); //kdotg_i;
|
||||
h_i=myH[ji]-myV[ji]*kk-outerProductSymm(kPoints[j],myG[jr]); //kdotg_r;
|
||||
for(int t=0; t<D*D; ++t)
|
||||
grad_grad_psi[psiIndex](t)= std::complex<TT>(c*h_r(t)-s*h_i(t), c*h_i(t)+s*h_r(t));
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void evaluate_vgh(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgh(r,0,psi,dpsi,grad_grad_psi);
|
||||
}
|
||||
|
||||
template<typename VGL>
|
||||
void evaluate_vgl_combo(const ParticleSet& P, const int iat, VGL& dpsi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
}
|
||||
};
|
||||
|
||||
/** adoptor class to match std::complex<ST> spline with TT real SPOs
|
||||
* @tparam ST precision of spline
|
||||
* @tparam TT precision of SPOs
|
||||
* @tparam D dimension
|
||||
*
|
||||
* Requires temporage storage and multiplication of phase vectors
|
||||
*/
|
||||
template<typename ST, typename TT, unsigned D>
|
||||
struct SplineC2RPackedAdoptor: public SplineAdoptorBase<ST,D>
|
||||
{
|
||||
typedef typename einspline_traits<ST,D>::SplineType SplineType;
|
||||
typedef typename einspline_traits<ST,D>::BCType BCType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::PointType PointType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::SingleSplineType SingleSplineType;
|
||||
|
||||
using SplineAdoptorBase<ST,D>::first_spo;
|
||||
using SplineAdoptorBase<ST,D>::last_spo;
|
||||
using SplineAdoptorBase<ST,D>::GGt;
|
||||
using SplineAdoptorBase<ST,D>::PrimLattice;
|
||||
using SplineAdoptorBase<ST,D>::kPoints;
|
||||
using SplineAdoptorBase<ST,D>::MakeTwoCopies;
|
||||
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myV;
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myL;
|
||||
typename OrbitalSetTraits<ST>::GradVector_t myG;
|
||||
typename OrbitalSetTraits<ST>::HessVector_t myH;
|
||||
typename OrbitalSetTraits<ST>::GradHessVector_t myGH;
|
||||
|
||||
///Actual spline table, multi_bspline_3d_(d,s)
|
||||
SplineType *MultiSpline;
|
||||
///number of points of the original grid
|
||||
int BaseN[3];
|
||||
///offset of the original grid, always 0
|
||||
int BaseOffset[3];
|
||||
|
||||
std::vector<ST> KdotR;
|
||||
std::vector<ST> CosV;
|
||||
std::vector<ST> SinV;
|
||||
std::vector<ST> mKK;
|
||||
std::vector<Tensor<ST,D> > KK; //k^k
|
||||
|
||||
SplineC2RPackedAdoptor():MultiSpline(0)
|
||||
{
|
||||
this->is_complex=true;
|
||||
this->AdoptorName="SplineC2RPackedAdoptor";
|
||||
this->KeyWord="C2RPacked";
|
||||
}
|
||||
|
||||
virtual ~SplineC2RPackedAdoptor() {}
|
||||
|
||||
inline void resizeStorage(int n, int nvals)
|
||||
{
|
||||
SplineAdoptorBase<ST,D>::init_base(n);
|
||||
myV.resize(2*n);
|
||||
myL.resize(2*n);
|
||||
myG.resize(2*n);
|
||||
myH.resize(2*n);
|
||||
CosV.resize(n);
|
||||
SinV.resize(n);
|
||||
KdotR.resize(n);
|
||||
}
|
||||
|
||||
template<typename GT, typename BCT>
|
||||
void create_spline(GT& xyz_g, BCT& xyz_bc)
|
||||
{
|
||||
resize_kk();
|
||||
MultiSpline=einspline::create(MultiSpline,xyz_g,xyz_bc,myV.size());
|
||||
for(int i=0; i<D; ++i)
|
||||
{
|
||||
BaseOffset[i]=0;
|
||||
BaseN[i]=xyz_g[i].num+3;
|
||||
}
|
||||
qmc_common.memory_allocated += MultiSpline->coefs_size*sizeof(ST);
|
||||
}
|
||||
|
||||
inline void resize_kk()
|
||||
{
|
||||
mKK.resize(kPoints.size());
|
||||
for(int i=0; i<kPoints.size(); ++i)
|
||||
mKK[i]=-dot(kPoints[i],kPoints[i]);
|
||||
KK.resize(kPoints.size());
|
||||
for(int i=0; i<kPoints.size(); ++i)
|
||||
KK[i]=outerProduct(kPoints[i],kPoints[i]);
|
||||
}
|
||||
|
||||
void set_spline(ST* restrict psi_r, ST* restrict psi_i, int twist, int ispline, int level)
|
||||
{
|
||||
//if(mKK.empty()) resize_kk();
|
||||
einspline::set(MultiSpline, 2*ispline, psi_r);
|
||||
einspline::set(MultiSpline, 2*ispline+1, psi_i);
|
||||
}
|
||||
|
||||
inline void set_spline(SingleSplineType* spline_r, SingleSplineType* spline_i, int twist, int ispline, int level)
|
||||
{
|
||||
//if(mKK.empty()) resize_kk();
|
||||
einspline::set(MultiSpline, 2*ispline,spline_r, BaseOffset, BaseN);
|
||||
einspline::set(MultiSpline, 2*ispline+1,spline_i, BaseOffset, BaseN);
|
||||
}
|
||||
|
||||
inline void set_spline_domain(SingleSplineType* spline_r, SingleSplineType* spline_i,
|
||||
int twist, int ispline, const int* offset_l, const int* mesh_l)
|
||||
{
|
||||
if(mKK.empty()) resize_kk();
|
||||
einspline::set(MultiSpline, 2*ispline, spline_r, offset_l, mesh_l);
|
||||
einspline::set(MultiSpline, 2*ispline+1,spline_i, offset_l, mesh_l);
|
||||
}
|
||||
|
||||
bool read_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.read(bigtable,o.str().c_str());//"spline_0");
|
||||
}
|
||||
|
||||
bool write_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.write(bigtable,o.str().c_str());//"spline_0");
|
||||
}
|
||||
|
||||
template<typename VV>
|
||||
inline void assign_v(const PointType& r, int bc_sign, VV& psi)
|
||||
{
|
||||
const int N=kPoints.size();
|
||||
for(int j=0; j<N; ++j)
|
||||
KdotR[j]=-dot(r,kPoints[j]);
|
||||
eval_e2iphi(N,KdotR.data(),CosV.data(),SinV.data());
|
||||
//for(int j=0; j<N; ++j) sincos(-dot(r,kPoints[j]),&SinV[j],&CosV[j]);
|
||||
int psiIndex=first_spo;
|
||||
for (int j=0,jr=0; j<N; j++,jr+=2)
|
||||
{
|
||||
psi[psiIndex] = static_cast<TT>(myV[jr]*CosV[j]-myV[jr+1]*SinV[j]);
|
||||
psiIndex++;
|
||||
if (MakeTwoCopies[j])
|
||||
{
|
||||
psi[psiIndex] = static_cast<TT>(myV[jr+1]*CosV[j]+myV[jr]*SinV[j]);
|
||||
psiIndex++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV>
|
||||
inline void evaluate_v(const ParticleSet& P, const int iat, VV& psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate(MultiSpline,ru,myV);
|
||||
assign_v(r,0,psi);
|
||||
}
|
||||
|
||||
template<typename VV, typename GV>
|
||||
inline void assign_vgl(const PointType& r, int bc_sign, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
const int N=kPoints.size();
|
||||
for (int j=0; j<2*N; j++)
|
||||
myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
for (int j=0; j<2*N; j++)
|
||||
myL[j] = trace(myH[j],GGt);
|
||||
const ST zero(0);
|
||||
const ST two(2);
|
||||
for(int j=0; j<N; ++j)
|
||||
KdotR[j]=-dot(r,kPoints[j]);
|
||||
eval_e2iphi(N,KdotR.data(),CosV.data(),SinV.data());
|
||||
//for(int j=0; j<N; ++j) sincos(-dot(r,kPoints[j]),&SinV[j],&CosV[j]);
|
||||
int psiIndex=first_spo;
|
||||
TinyVector<ST,D> g_r, g_i;
|
||||
for (int j=0,jr=0,ji=1; j<N; j++,jr+=2,ji+=2)
|
||||
{
|
||||
g_r=myG[jr]+myV[ji]*kPoints[j]; // \f$\nabla \psi_r + {\bf k}\psi_i\f$
|
||||
g_i=myG[ji]-myV[jr]*kPoints[j]; // \f$\nabla \psi_i - {\bf k}\psi_r\f$
|
||||
myL[jr]+=mKK[j]*myV[jr]+two*dot(kPoints[j],myG[ji]);
|
||||
myL[ji]+=mKK[j]*myV[ji]-two*dot(kPoints[j],myG[jr]);
|
||||
psi[psiIndex]=CosV[j]*myV[jr]-SinV[j]*myV[ji];
|
||||
dpsi[psiIndex]=CosV[j]*g_r-SinV[j]*g_i; // multiply phase
|
||||
d2psi[psiIndex]=CosV[j]*myL[jr]-SinV[j]*myL[ji];
|
||||
++psiIndex;
|
||||
if(MakeTwoCopies[j])
|
||||
{
|
||||
psi[psiIndex]=CosV[j]*myV[ji]+SinV[j]*myV[jr];
|
||||
dpsi[psiIndex]=CosV[j]*g_i+SinV[j]*g_r;
|
||||
d2psi[psiIndex]=CosV[j]*myL[ji]+SinV[j]*myL[jr];
|
||||
++psiIndex;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV, typename GV>
|
||||
inline void evaluate_vgl(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgl(r,0,psi,dpsi,d2psi);
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void assign_vgh(const PointType& r, int bc_sign, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
//convert to Cartesian G and Hess
|
||||
const int N=kPoints.size();
|
||||
for (int j=0; j<2*N; j++)
|
||||
myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
for (int j=0; j<2*N; j++)
|
||||
myH[j] = dot(myH[j],GGt);
|
||||
ST s,c;
|
||||
PointType g_r, g_i;
|
||||
Tensor<ST,D> h_r,h_i;
|
||||
//can easily make three independent loops
|
||||
int psiIndex=first_spo;
|
||||
for (int j=0,jr=0,ji=1; j<N; j++,jr+=2,ji+=2)
|
||||
{
|
||||
g_r=myG[jr]+myV[ji]*kPoints[j]; // \f$\nabla \psi_r + {\bf k}\psi_i\f$
|
||||
g_i=myG[ji]-myV[jr]*kPoints[j]; // \f$\nabla \psi_i - {\bf k}\psi_r\f$
|
||||
//kk=outerProduct(kPoints[j],kPoints[j]); // \f$kk=k^k \f$
|
||||
h_r=myH[jr]-myV[jr]*KK[j]+outerProductSymm(kPoints[j],myG[ji]);
|
||||
h_i=myH[ji]-myV[ji]*KK[j]-outerProductSymm(kPoints[j],myG[jr]);
|
||||
sincos(-dot(r,kPoints[j]),&s,&c); //e-ikr (beware of -1)
|
||||
psi[psiIndex]=c*myV[jr]-s*myV[ji];
|
||||
for(int idim=0; idim<D; ++idim)
|
||||
dpsi[psiIndex][idim]=c*g_r[idim]-s*g_i[idim];
|
||||
for(int t=0; t<D*D; ++t)
|
||||
grad_grad_psi[psiIndex](t)=c*h_r(t)-s*h_i(t);
|
||||
++psiIndex;
|
||||
if(MakeTwoCopies[j])
|
||||
{
|
||||
psi[psiIndex]=s*myV[jr]+c*myV[ji];
|
||||
for(int idim=0; idim<D; ++idim)
|
||||
dpsi[psiIndex][idim]=c*g_i[idim]+s*g_r[idim];
|
||||
for(int t=0; t<D*D; ++t)
|
||||
grad_grad_psi[psiIndex](t)=s*h_r(t)+c*h_i(t);
|
||||
++psiIndex;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void evaluate_vgh(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru(PrimLattice.toUnit_floor(r));
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgh(r,0,psi,dpsi,grad_grad_psi);
|
||||
}
|
||||
|
||||
template<typename VGL>
|
||||
void evaluate_vgl_combo(const ParticleSet& P, const int iat, VGL& dpsi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
// /** adoptor class to match std::complex<ST> spline with std::complex<TT> SPOs, just references
|
||||
// * @tparam ST precision of spline
|
||||
// * @tparam TT precision of SPOs
|
||||
// * @tparam D dimension
|
||||
// *
|
||||
// * This is a reference implementation but NOT USED. Use Packed version
|
||||
// */
|
||||
// template<typename ST, typename TT, unsigned D>
|
||||
// struct SplineC2CAdoptor
|
||||
// {
|
||||
// typedef ST real_type;
|
||||
// typedef std::complex<ST> value_type;
|
||||
// typedef typename einspline_traits<value_type,D>::SplineType SplineType;
|
||||
// typedef typename einspline_traits<value_type,D>::BCType BCType;
|
||||
// typedef typename OrbitalSetTraits<value_type>::ValueVector_t StorageValueVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::GradVector_t StorageGradVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::HessVector_t StorageHessVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::GradHessVector_t StorageGradHessVector_t;
|
||||
//
|
||||
// typedef CrystalLattice<ST,D> UnitCellType;
|
||||
// typedef TinyVector<ST,D> PointType;
|
||||
//
|
||||
// SplineType *MultiSpline;
|
||||
// UnitCellType SuperLattice;
|
||||
// UnitCellType PrimLattice;
|
||||
// std::vector<PointType> kPoints;
|
||||
// TinyVector<int,D> HalfG;
|
||||
// std::vector<bool> MakeTwoCopies;
|
||||
// Tensor<real_type,D> GGt;
|
||||
//
|
||||
// std::vector<real_type> phase;
|
||||
// std::vector<value_type> eikr;
|
||||
//
|
||||
// StorageValueVector_t myV;
|
||||
// StorageValueVector_t myL;
|
||||
// StorageGradVector_t myG;
|
||||
// StorageHessVector_t myH;
|
||||
// StorageGradHessVector_t myGH;
|
||||
//
|
||||
// inline void resizeStorage(int n, int nvals)
|
||||
// {
|
||||
// GGt=dot(transpose(PrimLattice.G),PrimLattice.G);
|
||||
//
|
||||
// kPoints.resize(n);
|
||||
// MakeTwoCopies.resize(n);
|
||||
// myV.resize(n);
|
||||
// myL.resize(n);
|
||||
// myG.resize(n);
|
||||
// myH.resize(n);
|
||||
// myH.resize(n);
|
||||
// }
|
||||
//
|
||||
// inline bool isready()
|
||||
// {
|
||||
// return true;
|
||||
// }
|
||||
//
|
||||
// template<typename VV>
|
||||
// inline void evaluate_v(const PointType& r, VV& psi)
|
||||
// {
|
||||
// PointType ru(PrimLattice.toUnit(r));
|
||||
// for (int i=0; i<D; i++) ru[i] -= std::floor (ru[i]);
|
||||
// einspline::evaluate(MultiSpline,ru,myV);
|
||||
//
|
||||
// const int N=myV.size();
|
||||
// //computePhases(r);
|
||||
// //simd::multiadd(psi.size(),eikr.data(),psi.data());
|
||||
// register ST s,c;
|
||||
// for(int j=0; j<psi.size(); ++j)
|
||||
// {
|
||||
// sincos(-dot(r,kPoints[j]),&s,&c);
|
||||
// psi[j]= std::complex<TT>(static_cast<TT>(c*myV[j].real()-s*myV[j].imag()),
|
||||
// static_cast<TT>(c*myV[j].imag()+s*myV[j].real()));
|
||||
// }
|
||||
//
|
||||
// }
|
||||
//
|
||||
// template<typename VV, typename GV>
|
||||
// inline void evaluate_vgl(const PointType& r, VV& psi, GV& dpsi, VV& d2psi)
|
||||
// {
|
||||
// PointType ru(PrimLattice.toUnit(r));
|
||||
// for (int i=0; i<D; i++) ru[i] -= std::floor (ru[i]);
|
||||
// einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
//
|
||||
// const int N=kPoints.size();
|
||||
// for (int j=0; j<N; j++) myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
// for (int j=0; j<N; j++) myL[j] = trace(myH[j],GGt);
|
||||
//
|
||||
// const ST two=2.0;
|
||||
// TinyVector<std::complex<ST>,D> ck;
|
||||
// register ST s,c;
|
||||
// for (int j=0; j<psi.size(); j++)
|
||||
// {
|
||||
// ST kk=dot(kPoints[j],kPoints[j]);
|
||||
// for(int i=0; i<D; ++i) ck[i]= std::complex<ST>(0.0,kPoints[j][i]);
|
||||
// sincos(-dot(r,kPoints[j]),&s,&c);
|
||||
// std::complex<ST> e_mikr(c,s);
|
||||
// convert(e_mikr * myV[j], psi[j]);
|
||||
// convert(e_mikr*(-myV[j]*ck + myG[j]), dpsi[j]);
|
||||
// convert(e_mikr*(-myV[j]*kk - two*dot(ck,myG[j]) + myL[j]), d2psi[j]);
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// template<typename VV, typename GV, typename GGV>
|
||||
// void evaluate_vgh(const PointType& r, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
// {
|
||||
// }
|
||||
// };
|
||||
//
|
||||
// /** adoptor class to match std::complex<ST> spline with TT real SPOs
|
||||
// * @tparam ST precision of spline
|
||||
// * @tparam TT precision of SPOs
|
||||
// * @tparam D dimension
|
||||
// *
|
||||
// * Requires temporage storage and multiplication of phase vectors
|
||||
// */
|
||||
// template<typename ST, typename TT, unsigned D>
|
||||
// struct SplineC2RAdoptor
|
||||
// {
|
||||
// typedef ST real_type;
|
||||
// typedef std::complex<ST> value_type;
|
||||
//
|
||||
// typedef typename einspline_traits<value_type,D>::SplineType SplineType;
|
||||
// typedef typename einspline_traits<value_type,D>::BCType BCType;
|
||||
// typedef typename OrbitalSetTraits<value_type>::ValueVector_t StorageValueVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::GradVector_t StorageGradVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::HessVector_t StorageHessVector_t;
|
||||
// typedef typename OrbitalSetTraits<value_type>::GradHessVector_t StorageGradHessVector_t;
|
||||
//
|
||||
// typedef CrystalLattice<ST,D> UnitCellType;
|
||||
// typedef TinyVector<ST,D> PointType;
|
||||
//
|
||||
// SplineType* MultiSpline;
|
||||
// UnitCellType SuperLattice;
|
||||
// UnitCellType PrimLattice;
|
||||
// TinyVector<int,D> HalfG;
|
||||
// Tensor<real_type,D> GGt;
|
||||
// std::vector<PointType> kPoints;
|
||||
// std::vector<bool> MakeTwoCopies;
|
||||
// std::vector<real_type> CosV;
|
||||
// std::vector<real_type> SinV;
|
||||
// std::vector<value_type> mKK;
|
||||
//
|
||||
// // Temporary storage for Eispline calls
|
||||
// StorageValueVector_t myV;
|
||||
// StorageValueVector_t myL;
|
||||
// StorageGradVector_t myG;
|
||||
// StorageHessVector_t myH;
|
||||
// StorageGradHessVector_t myGH;
|
||||
//
|
||||
// SplineC2RAdoptor():MultiSpline(0) { }
|
||||
//
|
||||
// virtual ~SplineC2RAdoptor(){}
|
||||
//
|
||||
// inline void resizeStorage(int n, int nvals)
|
||||
// {
|
||||
//
|
||||
// GGt=dot(transpose(PrimLattice.G),PrimLattice.G);
|
||||
//
|
||||
// HalfG=0;
|
||||
// kPoints.resize(n);
|
||||
// MakeTwoCopies.resize(n);
|
||||
//
|
||||
// myV.resize(n);
|
||||
// myL.resize(n);
|
||||
// myG.resize(n);
|
||||
// myH.resize(n);
|
||||
// CosV.resize(n);
|
||||
// SinV.resize(n);
|
||||
// }
|
||||
//
|
||||
// inline bool isready()
|
||||
// {
|
||||
// mKK.resize(kPoints.size());
|
||||
// for(int i=0; i<kPoints.size(); ++i) mKK[i]=-dot(kPoints[i],kPoints[i]);
|
||||
// return true;
|
||||
// }
|
||||
//
|
||||
// template<typename VV>
|
||||
// inline void evaluate_v(const PointType& r, VV& psi)
|
||||
// {
|
||||
// PointType ru(PrimLattice.toUnit(r));
|
||||
// for (int i=0; i<D; i++) ru[i] -= std::floor (ru[i]);
|
||||
// einspline::evaluate(MultiSpline,ru,myV);
|
||||
// const int N=kPoints.size();
|
||||
// for(int j=0; j<N; ++j) sincos(-dot(r,kPoints[j]),&SinV[j],&CosV[j]);
|
||||
//
|
||||
// int psiIndex = 0;
|
||||
// for (int j=0; j<N; j++)
|
||||
// {
|
||||
// psi[psiIndex] = static_cast<TT>(myV[j].real()*CosV[j]-myV[j].imag()*SinV[j]);
|
||||
// //psi[psiIndex] = static_cast<TT>(myV[j].real());
|
||||
// psiIndex++;
|
||||
// if (MakeTwoCopies[j])
|
||||
// {
|
||||
// //psi[psiIndex] = static_cast<TT>(myV[j].imag());
|
||||
// psi[psiIndex] = static_cast<TT>(myV[j].imag()*CosV[j]+myV[j].real()*SinV[j]);
|
||||
// psiIndex++;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// template<typename VV, typename GV>
|
||||
// inline void evaluate_vgl(const PointType& r, VV& psi, GV& dpsi, VV& d2psi)
|
||||
// {
|
||||
// PointType ru(PrimLattice.toUnit(r));
|
||||
// for (int i=0; i<D; i++) ru[i] -= std::floor (ru[i]);
|
||||
//
|
||||
// einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
// const int N=kPoints.size();
|
||||
// for (int j=0; j<N; j++) myG[j] = dot(PrimLattice.G, myG[j]);
|
||||
// for (int j=0; j<N; j++) myL[j] = trace(myH[j],GGt);
|
||||
//
|
||||
// const ST zero=0.0;
|
||||
// const ST two=2.0;
|
||||
// TinyVector<std::complex<ST>,D> ck;
|
||||
// ST s,c;
|
||||
// for(int j=0; j<N; ++j)
|
||||
// {
|
||||
// for (int n=0; n<D; n++)ck[n] = std::complex<ST>(zero,-kPoints[j][n]);
|
||||
// sincos (-dot(r,kPoints[j]), &s, &c);
|
||||
// std::complex<ST> e_mikr(c,s);
|
||||
// myV[j] = e_mikr*myV[j];
|
||||
// myL[j] = -dot(kPoints[j],kPoints[j])*myV[j]+e_mikr*(two*dot(ck,myG[j]) + myL[j]);
|
||||
// myG[j] = myV[j]*ck+e_mikr*myG[j];
|
||||
// }
|
||||
//
|
||||
// //const std::complex<ST> eye (0.0, 1.0);
|
||||
// //ST s,c;
|
||||
// //for(int j=0; j<N; ++j)
|
||||
// //{
|
||||
// // std::complex<ST> u = myV[j];
|
||||
// // TinyVector<std::complex<ST>,D> gradu = myG[j];
|
||||
// // std::complex<ST> laplu = myL[j];
|
||||
// // PointType k = kPoints[j];
|
||||
// // TinyVector<std::complex<ST>,D> ck;
|
||||
// // for (int n=0; n<D; n++) ck[n] = k[n];
|
||||
// // sincos (-dot(r,k), &s, &c);
|
||||
// // std::complex<ST> e_mikr (c,s);
|
||||
// // myV[j] = e_mikr*u;
|
||||
// // myG[j] = e_mikr*(-eye*u*ck + gradu);
|
||||
// // myL[j] = e_mikr*(-dot(k,k)*u - two*eye*dot(ck,gradu) + laplu);
|
||||
// //}
|
||||
//
|
||||
// int psiIndex=0;
|
||||
// for(int j=0; j<N; ++j)
|
||||
// {
|
||||
// psi[psiIndex]=static_cast<TT>(myV[j].real());
|
||||
// for(int n=0; n<D; ++n) dpsi[psiIndex][n]=static_cast<TT>(myG[j][n].real());
|
||||
// d2psi[psiIndex]=static_cast<TT>(myL[j].real());
|
||||
// ++psiIndex;
|
||||
// if(MakeTwoCopies[j])
|
||||
// {
|
||||
// psi[psiIndex]=static_cast<TT>(myV[j].imag());
|
||||
// for(int n=0; n<D; ++n) dpsi[psiIndex][n]=static_cast<TT>(myG[j][n].imag());
|
||||
// d2psi[psiIndex]=static_cast<TT>(myL[j].imag());
|
||||
// ++psiIndex;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// template<typename VV, typename GV, typename GGV>
|
||||
// void evaluate_vgh(const PointType& r, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
// {
|
||||
// }
|
||||
// };
|
||||
//
|
||||
}
|
||||
#endif
|
||||
|
|
@ -18,8 +18,8 @@
|
|||
*
|
||||
* The most general reader class for SplineAdoptor using the full single grid for the supercell
|
||||
* - SplineR2RAdoptor
|
||||
* - SplineC2CPackedAdoptor
|
||||
* - SplineC2RPackedAdoptor
|
||||
* - SplineC2CAdoptor
|
||||
* - SplineC2RAdoptor
|
||||
* Each band is initialized with UBspline_3d_d and both real and imaginary parts are passed to the adoptor object
|
||||
* which will convert the data type to their internal precision.
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -1,264 +0,0 @@
|
|||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
// This file is distributed under the University of Illinois/NCSA Open Source License.
|
||||
// See LICENSE file in top directory for details.
|
||||
//
|
||||
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
|
||||
//
|
||||
// File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_EINSPLINE_R2RADOPTOR_H
|
||||
#define QMCPLUSPLUS_EINSPLINE_R2RADOPTOR_H
|
||||
|
||||
#include <OhmmsSoA/Container.h>
|
||||
#include <spline2/MultiBspline.hpp>
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** adoptor class to match ST real spline with TT real SPOs
|
||||
* @tparam ST precision of spline
|
||||
* @tparam TT precision of SPOs
|
||||
* @tparam D dimension
|
||||
*/
|
||||
template<typename ST, typename TT, unsigned D>
|
||||
struct SplineR2RAdoptor: public SplineAdoptorBase<ST,D>
|
||||
{
|
||||
typedef typename einspline_traits<ST,D>::SplineType SplineType;
|
||||
typedef typename einspline_traits<ST,D>::BCType BCType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::PointType PointType;
|
||||
typedef typename SplineAdoptorBase<ST,D>::SingleSplineType SingleSplineType;
|
||||
|
||||
using SplineAdoptorBase<ST,D>::first_spo;
|
||||
using SplineAdoptorBase<ST,D>::last_spo;
|
||||
using SplineAdoptorBase<ST,D>::HalfG;
|
||||
using SplineAdoptorBase<ST,D>::GGt;
|
||||
using SplineAdoptorBase<ST,D>::PrimLattice;
|
||||
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myV;
|
||||
typename OrbitalSetTraits<ST>::ValueVector_t myL;
|
||||
typename OrbitalSetTraits<ST>::GradVector_t myG;
|
||||
typename OrbitalSetTraits<ST>::HessVector_t myH;
|
||||
typename OrbitalSetTraits<ST>::GradHessVector_t myGH;
|
||||
|
||||
SplineType *MultiSpline;
|
||||
|
||||
///number of points of the original grid
|
||||
int BaseN[3];
|
||||
///offset of the original grid, always 0
|
||||
int BaseOffset[3];
|
||||
|
||||
SplineR2RAdoptor(): MultiSpline(0)
|
||||
{
|
||||
this->is_complex=false;
|
||||
this->AdoptorName="SplineR2RAdoptor";
|
||||
this->KeyWord="R2R";
|
||||
}
|
||||
|
||||
void resizeStorage(int n, int nv)
|
||||
{
|
||||
SplineAdoptorBase<ST,D>::init_base(n);
|
||||
myV.resize(n);
|
||||
myL.resize(n);
|
||||
myG.resize(n);
|
||||
myH.resize(n);
|
||||
myGH.resize(n);
|
||||
}
|
||||
|
||||
template<typename GT, typename BCT>
|
||||
void create_spline(GT& xyz_g, BCT& xyz_bc)
|
||||
{
|
||||
GGt=dot(transpose(PrimLattice.G),PrimLattice.G);
|
||||
MultiSpline=einspline::create(MultiSpline,xyz_g,xyz_bc,myV.size());
|
||||
|
||||
qmc_common.memory_allocated += MultiSpline->coefs_size*sizeof(ST);
|
||||
|
||||
for(int i=0; i<D; ++i)
|
||||
{
|
||||
BaseOffset[i]=0;
|
||||
BaseN[i]=xyz_g[i].num+3;
|
||||
}
|
||||
}
|
||||
|
||||
void create_spline(TinyVector<int,D>& mesh, int n)
|
||||
{
|
||||
Ugrid xyz_grid[D];
|
||||
BCType xyz_bc[D];
|
||||
for(int i=0; i<D; ++i)
|
||||
{
|
||||
xyz_grid[i].start = 0.0;
|
||||
xyz_grid[i].end = 1.0;
|
||||
xyz_grid[i].num = mesh[i];
|
||||
xyz_bc[i].lCode=xyz_bc[i].rCode=(HalfG[i])? ANTIPERIODIC:PERIODIC;
|
||||
BaseOffset[i]=0;
|
||||
BaseN[i]=xyz_grid[i].num+3;
|
||||
}
|
||||
SplineType* dummy=0;
|
||||
MultiSpline=einspline::create(dummy,xyz_grid,xyz_bc,n);
|
||||
|
||||
qmc_common.memory_allocated += MultiSpline->coefs_size*sizeof(ST);
|
||||
}
|
||||
|
||||
inline bool isready()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
inline void set_spline(ST* restrict psi_r, ST* restrict psi_i, int twist, int ispline, int level)
|
||||
{
|
||||
einspline::set(MultiSpline, ispline,psi_r);
|
||||
}
|
||||
|
||||
inline void set_spline(SingleSplineType* spline_r, SingleSplineType* spline_i, int twist, int ispline, int level)
|
||||
{
|
||||
einspline::set(MultiSpline, ispline,spline_r, BaseOffset,BaseN);
|
||||
}
|
||||
|
||||
inline void set_spline_domain(SingleSplineType* spline_r, SingleSplineType* spline_i,
|
||||
int twist, int ispline, const int* offset_l, const int* mesh_l)
|
||||
{
|
||||
einspline::set(MultiSpline, ispline, spline_r, offset_l, mesh_l);
|
||||
}
|
||||
|
||||
bool read_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.read(bigtable,o.str().c_str()); //"spline_0");
|
||||
}
|
||||
|
||||
bool write_splines(hdf_archive& h5f)
|
||||
{
|
||||
std::ostringstream o;
|
||||
o<<"spline_" << SplineAdoptorBase<ST,D>::MyIndex;
|
||||
einspline_engine<SplineType> bigtable(MultiSpline);
|
||||
return h5f.write(bigtable,o.str().c_str());//"spline_0");
|
||||
}
|
||||
|
||||
/** convert postion in PrimLattice unit and return sign */
|
||||
inline int convertPos(const PointType& r, PointType& ru)
|
||||
{
|
||||
ru=PrimLattice.toUnit(r);
|
||||
int bc_sign=0;
|
||||
for(int i=0; i<D; i++)
|
||||
if( -std::numeric_limits<ST>::epsilon() < ru[i] && ru[i] < 0 )
|
||||
ru[i] = ST(0.0);
|
||||
else
|
||||
{
|
||||
ST img = std::floor(ru[i]);
|
||||
ru[i] -= img;
|
||||
bc_sign += HalfG[i] * (int)img;
|
||||
}
|
||||
return bc_sign;
|
||||
}
|
||||
|
||||
/** assign myV to psi
|
||||
*/
|
||||
template<typename VV>
|
||||
inline void assign_v(const PointType& r, int bc_sign, VV& psi)
|
||||
{
|
||||
if (bc_sign & 1)
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=static_cast<TT>(-myV[j]);
|
||||
else
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=static_cast<TT>(myV[j]);
|
||||
}
|
||||
|
||||
template<typename VV>
|
||||
inline void evaluate_v(const ParticleSet& P, const int iat, VV& psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru;
|
||||
int bc_sign=convertPos(r,ru);
|
||||
einspline::evaluate(MultiSpline,ru,myV);
|
||||
assign_v(r,bc_sign,psi);
|
||||
}
|
||||
|
||||
/** assign internal data to psi's
|
||||
*/
|
||||
template<typename VV, typename GV>
|
||||
inline void assign_vgl(const PointType& r, int bc_sign, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
if (bc_sign & 1)
|
||||
{
|
||||
const ST minus_one=-1.0;
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=-myV[j];
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
dpsi[psiIndex]=minus_one*dot(PrimLattice.G,myG[j]);
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
d2psi[psiIndex]=-trace(myH[j],GGt);
|
||||
}
|
||||
else
|
||||
{
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=myV[j];
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
dpsi[psiIndex]=dot(PrimLattice.G,myG[j]);
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
d2psi[psiIndex]=trace(myH[j],GGt);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV, typename GV>
|
||||
inline void evaluate_vgl(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, VV& d2psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru;
|
||||
int bc_sign=convertPos(r,ru);
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgl(r,bc_sign,psi,dpsi,d2psi);
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void assign_vgh(const PointType& r, int bc_sign, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
if (bc_sign & 1)
|
||||
{
|
||||
const ST minus_one=-1.0;
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=-myV[j];
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
dpsi[psiIndex]=minus_one*dot(PrimLattice.G,myG[j]);
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
grad_grad_psi[psiIndex]=minus_one*dot(myH[j],GGt);
|
||||
}
|
||||
else
|
||||
{
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
psi[psiIndex]=myV[j];
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
dpsi[psiIndex]=dot(PrimLattice.G,myG[j]);
|
||||
for(int psiIndex=first_spo,j=0; psiIndex<last_spo; ++psiIndex,++j)
|
||||
grad_grad_psi[psiIndex]=dot(myH[j],GGt);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename VV, typename GV, typename GGV>
|
||||
void evaluate_vgh(const ParticleSet& P, const int iat, VV& psi, GV& dpsi, GGV& grad_grad_psi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
PointType ru;
|
||||
int bc_sign=convertPos(r,ru);
|
||||
einspline::evaluate_vgh(MultiSpline,ru,myV,myG,myH);
|
||||
assign_vgh(r,bc_sign,psi,dpsi,grad_grad_psi);
|
||||
}
|
||||
|
||||
template<typename VGL>
|
||||
void evaluate_vgl_combo(const ParticleSet& P, const int iat, VGL& dpsi)
|
||||
{
|
||||
const PointType& r=P.R[iat];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -132,6 +132,47 @@ const char *particles =
|
|||
SPOSetBase *spo = einSet.createSPOSetFromXML(ein1);
|
||||
REQUIRE(spo != NULL);
|
||||
|
||||
#if !defined(QMC_CUDA) || defined(QMC_COMPLEX)
|
||||
// due to the different ordering of bands skip the tests on CUDA+Real builds
|
||||
// checking evaluations, reference values are not independently generated.
|
||||
// for vgl
|
||||
SPOSetBase::ValueMatrix_t psiM(elec_.R.size(),spo->getOrbitalSetSize());
|
||||
SPOSetBase::GradMatrix_t dpsiM(elec_.R.size(),spo->getOrbitalSetSize());
|
||||
SPOSetBase::ValueMatrix_t d2psiM(elec_.R.size(),spo->getOrbitalSetSize());
|
||||
spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM, dpsiM, d2psiM);
|
||||
|
||||
// value
|
||||
REQUIRE(psiM[1][0] == ComplexApprox(-0.8886948824).compare_real_only());
|
||||
REQUIRE(psiM[1][1] == ComplexApprox(1.4194120169).compare_real_only());
|
||||
// grad
|
||||
REQUIRE(dpsiM[1][0][0] == ComplexApprox(-0.0000183403).compare_real_only());
|
||||
REQUIRE(dpsiM[1][0][1] == ComplexApprox(0.1655139178).compare_real_only());
|
||||
REQUIRE(dpsiM[1][0][2] == ComplexApprox(-0.0000193077).compare_real_only());
|
||||
REQUIRE(dpsiM[1][1][0] == ComplexApprox(-1.3131694794).compare_real_only());
|
||||
REQUIRE(dpsiM[1][1][1] == ComplexApprox(-1.1174004078).compare_real_only());
|
||||
REQUIRE(dpsiM[1][1][2] == ComplexApprox(-0.8462534547).compare_real_only());
|
||||
// lapl
|
||||
REQUIRE(d2psiM[1][0] == ComplexApprox(1.3313053846).compare_real_only());
|
||||
REQUIRE(d2psiM[1][1] == ComplexApprox(-4.712583065).compare_real_only());
|
||||
|
||||
// for vgh
|
||||
SPOSetBase::ValueVector_t psiV(psiM[1],spo->getOrbitalSetSize());
|
||||
SPOSetBase::GradVector_t dpsiV(dpsiM[1],spo->getOrbitalSetSize());
|
||||
SPOSetBase::HessVector_t ddpsiV(spo->getOrbitalSetSize());
|
||||
spo->evaluate(elec_, 1, psiV, dpsiV, ddpsiV);
|
||||
|
||||
//hess
|
||||
REQUIRE(ddpsiV[1](0,0) == ComplexApprox(-2.3160984034).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](0,1) == ComplexApprox(1.8089479397).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](0,2) == ComplexApprox(0.5608575749).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](1,0) == ComplexApprox(1.8089479397).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](1,1) == ComplexApprox(-0.0799493616).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](1,2) == ComplexApprox(0.5237969314).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](2,0) == ComplexApprox(0.5608575749).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](2,1) == ComplexApprox(0.5237969314).compare_real_only());
|
||||
REQUIRE(ddpsiV[1](2,2) == ComplexApprox(-2.316497764).compare_real_only());
|
||||
#endif
|
||||
|
||||
#if 0
|
||||
// Dump values of the orbitals
|
||||
int orbSize= spo->getOrbitalSetSize();
|
||||
|
|
|
|||
|
|
@ -62,6 +62,14 @@ namespace qmcplusplus
|
|||
return h00*gg[0]+h01*gg[1]+h02*gg[2]+h11*gg[3]+h12*gg[4]+h22*gg[5];
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T v_m_v(T h00, T h01, T h02, T h11, T h12, T h22, T g1x, T g1y, T g1z, T g2x, T g2y, T g2z)
|
||||
{
|
||||
return g1x*g2x*h00+g1x*g2y*h01+g1x*g2z*h02
|
||||
+g1y*g2x*h01+g1y*g2y*h11+g1y*g2z*h12
|
||||
+g1z*g2x*h02+g1z*g2y*h12+g1z*g2z*h22;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
struct MultiBsplineData
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
|
||||
IF (NOT MIXED_PRECISION)
|
||||
SIMPLE_RUN_AND_CHECK(short-bccH_2x2x2_ae-grad-lap
|
||||
SIMPLE_RUN_AND_CHECK(short-bccH_2x2x2_ae-grad_lap
|
||||
"${CMAKE_SOURCE_DIR}/tests/solids/bccH_2x2x2_ae"
|
||||
grad_lap.xml
|
||||
1 1
|
||||
|
|
|
|||
|
|
@ -50,7 +50,7 @@
|
|||
</group>
|
||||
</particleset>
|
||||
<wavefunction name="psi0" target="e">
|
||||
<sposet_builder type="bspline" href="pwscf.pwscf.h5" tilematrix="2 0 0 0 2 0 0 0 2" twistnum="0" source="ion0" version="0.10" meshfactor="1.0" use_old_spline="yes">
|
||||
<sposet_builder type="bspline" href="pwscf.pwscf.h5" tilematrix="2 0 0 0 2 0 0 0 2" twistnum="0" source="ion0" version="0.10" meshfactor="1.0">
|
||||
<sposet type="bspline" name="spo_ud" size="8" spindataset="0"/>
|
||||
</sposet_builder>
|
||||
<determinantset>
|
||||
|
|
|
|||
|
|
@ -90,7 +90,7 @@
|
|||
</group>
|
||||
</particleset>
|
||||
<wavefunction name="psi0" target="e">
|
||||
<sposet_builder type="bspline" href="pwscf.pwscf.h5" tilematrix="3 0 0 0 3 0 0 0 3" twistnum="0" source="ion0" size="27" use_old_spline="yes">
|
||||
<sposet_builder type="bspline" href="pwscf.pwscf.h5" tilematrix="3 0 0 0 3 0 0 0 3" twistnum="0" source="ion0" size="27">
|
||||
<sposet type="bspline" name="spo_ud" size="27" spindataset="0"/>
|
||||
</sposet_builder>
|
||||
<determinantset>
|
||||
|
|
|
|||
Loading…
Reference in New Issue