scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Real orbitals now make use of +k/-k pairs. 

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@2565 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Ken Esler 2008-04-04 22:31:21 +00:00
parent 9be2456180
commit e0936987ac
4 changed files with 100 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
src/QMCWaveFunctions/EinsplineSet.cpp
View File

@ -288,11 +288,6 @@ namespace qmcplusplus {
EinsplineSetExtended<StorageType>::setOrbitalSetSize(int norbs)
{
OrbitalSetSize = norbs;
eikr.resize(norbs);
phase.resize(norbs);
StorageValueVector.resize(norbs);
StorageGradVector.resize(norbs);
StorageHessVector.resize(norbs);
}
template<typename StorageType> void
@ -304,14 +299,23 @@ namespace qmcplusplus {
for (int i=0; i<OHMMS_DIM; i++)
ru[i] -= std::floor (ru[i]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector);
computePhaseFactors(r);
for (int i=0; i<psi.size(); i++) {
//computePhaseFactors(r);
int N = StorageValueVector.size();
int psiIndex = 0;
for (int i=0; i<N; i++) {
PosType k = kPoints[i];
double s,c;
double phase = -dot(r, k);
sincos (phase, &s, &c);
complex<double> e_mikr (c,s);
convert (e_mikr*StorageValueVector[i], psi[i]);
complex<double> psi_val = e_mikr*StorageValueVector[i];
psi[psiIndex] = real(psi_val);
psiIndex++;
if (MakeTwoCopies[i]) {
psi[psiIndex] = imag(psi_val);
psiIndex++;
}
//convert (e_mikr*StorageValueVector[i], psi[i]);
}
}
@ -325,7 +329,7 @@ namespace qmcplusplus {
for (int i=0; i<OHMMS_DIM; i++)
ru[i] -= std::floor (ru[i]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector);
computePhaseFactors(r);
//computePhaseFactors(r);
for (int i=0; i<psi.size(); i++) {
PosType k = kPoints[i];
double s,c;
@ -362,7 +366,7 @@ namespace qmcplusplus {
ru[i] -= std::floor (ru[i]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
StorageGradVector, StorageHessVector);
computePhaseFactors(r);
//computePhaseFactors(r);
complex<double> eye (0.0, 1.0);
for (int j=0; j<psi.size(); j++) {
complex<double> u, laplu;
@ -397,7 +401,7 @@ namespace qmcplusplus {
ru[i] -= std::floor (ru[i]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
StorageGradVector, StorageHessVector);
computePhaseFactors(r);
//computePhaseFactors(r);
complex<double> eye (0.0, 1.0);
for (int j=0; j<psi.size(); j++) {
complex<double> u, laplu;
@ -451,9 +455,11 @@ namespace qmcplusplus {
ru[n] -= std::floor (ru[n]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
StorageGradVector, StorageHessVector);
computePhaseFactors(r);
//computePhaseFactors(r);
complex<double> eye (0.0, 1.0);
for (int j=0; j<OrbitalSetSize; j++) {
int N = StorageValueVector.size();
int psiIndex = 0;
for (int j=0; j<N; j++) {
complex<double> u, laplu;
TinyVector<complex<double>, OHMMS_DIM> gradu;
u = StorageValueVector[j];
@ -467,14 +473,35 @@ namespace qmcplusplus {
double phase = -dot(r, k);
sincos (phase, &s, &c);
complex<double> e_mikr (c,s);
psi(j,i) = real(e_mikr*u);
complex<double> psi_val = e_mikr*u;
TinyVector<complex<double>,OHMMS_DIM> psi_grad =
e_mikr*(-eye*u*ck + gradu);
complex<double> psi_lapl = e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu);
psi(psiIndex,i) = real(psi_val);
for (int n=0; n<3; n++)
dpsi(i,j) = real(e_mikr*(-eye*u*ck + gradu));
d2psi(i,j) = real(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu));
dpsi(i,psiIndex)[n] = real(psi_grad[n]);
d2psi(i,psiIndex) = real(psi_lapl);
psiIndex++;
if (MakeTwoCopies[j]) {
psi(psiIndex,i) = imag(psi_val);
for (int n=0; n<3; n++)
dpsi(i,psiIndex)[n] = imag(psi_grad[n]);
d2psi(i,psiIndex) = imag(psi_lapl);
psiIndex++;
}
// psi(j,i) = real(e_mikr*u);
// for (int n=0; n<3; n++)
// dpsi(i,j)[n] = real(e_mikr*(-eye*u*ck + gradu));
// d2psi(i,j) = real(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu));
//convert(e_mikr * u, psi(j,i));
//convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
//convert(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu), d2psi(i,j));
}
//cerr << "VGH psiIndex = " << psiIndex << endl;
}
}
@ -492,7 +519,7 @@ namespace qmcplusplus {
ru[n] -= std::floor (ru[n]);
EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
StorageGradVector, StorageHessVector);
computePhaseFactors(r);
//computePhaseFactors(r);
complex<double> eye (0.0, 1.0);
for (int j=0; j<OrbitalSetSize; j++) {
complex<double> u, laplu;

6
src/QMCWaveFunctions/EinsplineSet.h
View File

@ -179,7 +179,7 @@ namespace qmcplusplus {
StorageGradVector_t StorageGradVector;
StorageHessVector_t StorageHessVector;
// True if we should unpack this orbital into two copies
Vector<bool> MakeTwoCopies;
vector<bool> MakeTwoCopies;
// k-points for each orbital
Vector<TinyVector<double,OHMMS_DIM> > kPoints;
@ -221,12 +221,12 @@ namespace qmcplusplus {
inline void EinsplineSetExtended<StorageType>::computePhaseFactors(TinyVector<double,OHMMS_DIM> r)
{
#ifdef HAVE_MKL
for (int i=0; i<OrbitalSetSize; i++)
for (int i=0; i<kPoints.size(); i++)
phase[i] = -dot(r, kPoints[i]);
vzCIS(OrbitalSetSize, phase, (double*)eikr.data());
#else
double s, c;
for (int i=0; i<OrbitalSetSize; i++) {
for (int i=0; i<kPoints.size(); i++) {
phase[i] = -dot(r, kPoints[i]);
sincos (phase[i], &s, &c);
eikr[i] = complex<double>(c,s);

60
src/QMCWaveFunctions/EinsplineSetBuilder.cpp
View File

@ -464,8 +464,10 @@ namespace qmcplusplus {
for (int i=0; i<superSets[TwistNum].size(); i++)
IncludeTwists.push_back(superSets[TwistNum][i]);
vector<int> copyTwists;
// Now, find out which twists are distinct
DistinctTwists.clear();
#ifndef QMC_COMPLEX
vector<int> copyTwists;
for (int i=0; i<IncludeTwists.size(); i++) {
int ti = IncludeTwists[i];
PosType twist_i = TwistAngles[ti];
@ -475,11 +477,8 @@ namespace qmcplusplus {
PosType twist_j = TwistAngles[tj];
PosType sum = twist_i + twist_j;
PosType diff = twist_i - twist_j;
if (TwistPair (twist_i, twist_j)) {
cerr << "twist " << twist_i << " is paired with twist "
<< twist_j << endl;
if (TwistPair (twist_i, twist_j))
distinct = false;
}
}
if (distinct)
DistinctTwists.push_back(ti);
@ -501,6 +500,15 @@ namespace qmcplusplus {
fprintf (stderr, "Using %d copies of twist angle [%6.3f, %6.3f, %6.3f]\n",
MakeTwoCopies[i] ? 2 : 1, twist_i[0], twist_i[1], twist_i[2]);
}
#else
DistinctTwists.resize(IncludeTwists.size());
MakeTwoCopies.resize(IncludeTwists.size());
for (int i=0; i<IncludeTwists.size(); i++) {
DistinctTwists[i] = IncludeTwists[i];
MakeTwoCopies[i] = false;
}
#endif
}
// This function analyzes the twist vectors to see if they lay on a
@ -624,12 +632,13 @@ namespace qmcplusplus {
eigenstatesGroup = "/eigenstates";
SortBands.clear();
for (int ti=0; ti<IncludeTwists.size(); ti++) {
int tindex = IncludeTwists[ti];
for (int ti=0; ti<DistinctTwists.size(); ti++) {
int tindex = DistinctTwists[ti];
for (int bi=0; bi<NumBands; bi++) {
BandInfo band;
band.TwistIndex = tindex;
band.BandIndex = bi;
band.MakeTwoCopies = MakeTwoCopies[ti];
// Read eigenenergy from file
ostringstream ePath, sPath;
@ -664,6 +673,18 @@ namespace qmcplusplus {
cerr << "Sorting the bands now:\n";
sort (SortBands.begin(), SortBands.end());
}
int orbIndex = 0;
int numOrbs = 0;
while (numOrbs < OrbitalSet->getOrbitalSetSize()) {
if (SortBands[orbIndex].MakeTwoCopies)
numOrbs += 2;
else
numOrbs++;
orbIndex++;
}
NumDistinctOrbitals = orbIndex;
cerr << "We will read " << NumDistinctOrbitals << " distinct orbitals.\n";
}
@ -749,13 +770,21 @@ namespace qmcplusplus {
EinsplineSetBuilder::ReadBands
(int spin, EinsplineSetExtended<double>* orbitalSet)
{
int N = orbitalSet->getOrbitalSetSize();
//int N = orbitalSet->getOrbitalSetSize();
int N = NumDistinctOrbitals;
orbitalSet->kPoints.resize(N);
orbitalSet->MakeTwoCopies.resize(N);
orbitalSet->StorageValueVector.resize(N);
orbitalSet->StorageGradVector.resize(N);
orbitalSet->StorageHessVector.resize(N);
orbitalSet->phase.resize(N);
orbitalSet->eikr.resize(N);
// Read in k-points
orbitalSet->kPoints.resize(orbitalSet->getOrbitalSetSize());
for (int iorb=0; iorb<N; iorb++) {
int ti = SortBands[iorb].TwistIndex;
PosType twist = TwistAngles[ti];
orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
orbitalSet->MakeTwoCopies[iorb] = SortBands[iorb].MakeTwoCopies;
}
// First, check to see if we have already read this in
@ -876,13 +905,22 @@ namespace qmcplusplus {
EinsplineSetBuilder::ReadBands
(int spin, EinsplineSetExtended<complex<double> >* orbitalSet)
{
int N = orbitalSet->getOrbitalSetSize();
// int N = orbitalSet->getOrbitalSetSize();
int N = NumDistinctOrbitals;
cerr << "NumDistinctOrbitals = " << N << endl;
orbitalSet->kPoints.resize(N);
orbitalSet->MakeTwoCopies.resize(N);
orbitalSet->StorageValueVector.resize(N);
orbitalSet->StorageGradVector.resize(N);
orbitalSet->StorageHessVector.resize(N);
orbitalSet->phase.resize(N);
orbitalSet->eikr.resize(N);
// Read in k-points
orbitalSet->kPoints.resize(orbitalSet->getOrbitalSetSize());
for (int iorb=0; iorb<N; iorb++) {
int ti = SortBands[iorb].TwistIndex;
PosType twist = TwistAngles[ti];
orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
orbitalSet->MakeTwoCopies[iorb] = SortBands[iorb].MakeTwoCopies;
}
// First, check to see if we have already read this in

4
src/QMCWaveFunctions/EinsplineSetBuilder.h
View File

@ -80,6 +80,9 @@ namespace qmcplusplus {
struct BandInfo {
int TwistIndex, BandIndex, Spin;
double Energy;
// This is true if we should make distinct copies
// represeninting a +k, -k pair
bool MakeTwoCopies;
inline bool operator<(BandInfo other) const
{ return Energy < other.Energy; }
};
@ -142,6 +145,7 @@ namespace qmcplusplus {
// clone
vector<TinyVector<int,OHMMS_DIM> > UseTwists;
vector<int> IncludeTwists, DistinctTwists;
int NumDistinctOrbitals;
// This is true if the corresponding twist in DistinctTwists should
// should be used to generate two distinct orbitals from the real and
// imaginary parts.