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

several small changes in multidet implementation including: added ORMAS output to converter, fixed bug with 6tuple excitations, added CSF to the slow MSD class, etc Also added some funcitonality to cusp correction classes. 

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@5103 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Miguel Morales 2011-01-28 02:30:52 +00:00
parent d4e865921f
commit 6363406263
14 changed files with 678 additions and 491 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -20,7 +20,7 @@ SET(QMCPLUSPLUS_VERSION "${QMCPLUSPLUS_VERSION_MAJOR}.${QMCPLUSPLUS_VERSION_MINO
######################################################################
# Build level
######################################################################
SET(QMC_BUILD_LEVEL 1 CACHE INTEGER
SET(QMC_BUILD_LEVEL 3 CACHE INTEGER
"QMC Build Level: 1=bare, 2=developer, 3=experimental")
######################################################################
@ -145,7 +145,7 @@ CHECK_INCLUDE_FILE(unistd.h HAVE_UNISTD_H)
# MPIP_PROFILE profile mpi performance
######################################################################
SET(QMC_BUILD_STATIC 0 CACHE BOOL "Build static libraries and binary")
SET(BUILD_QMCTOOLS 0 CACHE BOOL "Build tools for QMCPACK")
SET(BUILD_QMCTOOLS 1 CACHE BOOL "Build tools for QMCPACK")
SET(BUILD_SANDBOX 0 CACHE BOOL "Build snadbox for testing")
SET(MPIP_PROFILE 0 CACHE BOOL "Build with mpip for mpi profile")

1
src/OhmmsPETE/OhmmsMatrix.h
View File

@ -36,6 +36,7 @@ public:
typedef const T* const_pointer;
typedef C Container_t;
typedef typename C::size_type size_type;
typedef typename Container_t::iterator iterator;
typedef Matrix<T,C> This_t;
Matrix():D1(0),D2(0){ } // Default Constructor initializes to zero.

70
src/QMCTools/GamesAsciiParser.cpp
View File

@ -175,8 +175,11 @@ void GamesAsciiParser::parse(const std::string& fname) {
if(multideterminant) {
fin.open(outputFile.c_str());
pivot_begin= fin.tellg();
//cout<<"looking for dets " <<endl;
//cout.flush();
if(lookFor(fin,"GUGA DISTINCT ROW TABLE")) {
cout<<"Found GUGA ROW TABLE, reading CSF." <<endl;
//cout.flush();
if(!lookFor(fin,"SYMMETRIES FOR THE",aline)) {
cerr<<"Could not find number of frozen core orbitals in output file.\n";
abort();
@ -187,19 +190,23 @@ void GamesAsciiParser::parse(const std::string& fname) {
NTOT=NEXT+NAC;
cout<<"# core, #active, #external: " <<NFZC <<" " <<NAC <<" " <<NEXT <<endl;
}
//cout.flush();
fin.seekg(pivot_begin);
getCSF(fin);
} else {
cout<<"Could not find GUGA ROW TABLE, looking for Slater Dets." <<endl;
//cout.flush();
fin.close(); fin.open(outputFile.c_str());
pivot_begin= fin.tellg();
if(lookFor(fin,"DIRECT DETERMINANT ORMAS-CI")) {
cout<<"Found ORMAS-CI" <<endl;
//cout.flush();
fin.close(); fin.open(outputFile.c_str());
pivot_begin= fin.tellg();
getORMAS(fin);
} else {
cout<<"Assuming ALDET-CI" <<endl;
//cout.flush();
fin.close(); fin.open(outputFile.c_str());
pivot_begin= fin.tellg();
getCI(fin);
@ -272,7 +279,8 @@ void GamesAsciiParser::getGeometry(std::istream& is) {
core.resize(NumberOfAtoms);
getwords(currentWords,is); // -------------
// this only works if all atoms have an ECP, fix later
for(int i=0; i<NumberOfAtoms; i++) {
// for(int i=0; i<NumberOfAtoms; i++) {
// fixing this problem
bool done=false;
while(!done) {
if(is.eof()) {
@ -282,45 +290,58 @@ void GamesAsciiParser::getGeometry(std::istream& is) {
getwords(currentWords,is);
if(currentWords.size() == 0) continue;
if(currentWords.size() >= 4) {
if(currentWords[0] == "THE" &&
currentWords[1] == "ECP" &&
currentWords[2] == "RUN" &&
currentWords[3] == "REMOVES" ) { done=true; }
}
if(currentWords[0] == "PARAMETERS" &&
currentWords[1] == "FOR" ) {
done=true;
std::vector<std::string>::iterator it;
//done=true;
std::vector<std::string>::iterator it,it0;
it = find(currentWords.begin(),currentWords.end(),"ZCORE");
it0 = find(currentWords.begin(),currentWords.end(),"ATOM");
if(it0 == currentWords.end()) {
cerr<<"Problem with ECP data. Didn't found ATOM tag\n";
cerr<<is <<endl;
abort();
}
it0++;
int nq0 = atoi(it0->c_str())-1;
if(it != currentWords.end()) {
it++;
core[i] = atoi(it->c_str());
q[i] -= core[i];
cout<<"Found ECP for atom " <<i <<" with zcore " <<core[i] <<endl;
core[nq0] = atoi(it->c_str());
q[nq0] -= core[nq0];
cout<<"Found ECP for atom " <<nq0 <<" with zcore " <<core[nq0] <<endl;
} else {
it = find(currentWords.begin(),currentWords.end(),"ATOM");
if(it == currentWords.end()) {
cerr<<"Problem with ECP data. Didn't found ATOM tag\n";
cerr<<"Atom: " <<i <<endl;
cerr<<"Atom: " <<nq0 <<endl;
abort();
}
std::vector<std::string>::iterator it2=it;
it2++;
int nq = atoi(it2->c_str());
if(nq != i+1) {
if(nq != nq0+1) {
cerr<<"Problem with ECP data. ID's don't agree\n";
cerr<<"Atom: " <<i <<endl;
cerr<<"Atom: " <<nq0 <<endl;
abort();
}
it = find(it2,currentWords.end(),"ATOM");
if(it == currentWords.end()) {
cerr<<"Problem with ECP data (2).\n";
cerr<<"Atom: " <<i <<endl;
cerr<<"Atom: " <<nq0 <<endl;
abort();
}
nq = atoi((it+1)->c_str());
core[i] = core[nq-1];
q[i] -= core[i];
cout<<"Found ECP for atom " <<i <<" with zcore " <<core[i] <<endl;
core[nq0] = core[nq-1];
q[nq0] -= core[nq0];
cout<<"Found ECP for atom " <<nq0 <<" with zcore " <<core[nq0] <<endl;
}
}
}
}
notfound=false;
} else {
if(currentWords.size() < 3 ) continue;
@ -761,8 +782,11 @@ void GamesAsciiParser::getCSF(std::istream& is)
int nq = atoi(currentWords[0].c_str());
pair<int,double> cic(nq,cof);
coeff2csf.push_back(cic);
if(NTOT <= 50) {
if(NTOT < 50) {
CSFocc.push_back(currentWords[2]);
} else if(NTOT == 50) {
CSFocc.push_back(currentWords[2]);
getwords(currentWords,is);
} else {
string tmp=currentWords[2];
getwords(currentWords,is);
@ -771,7 +795,7 @@ void GamesAsciiParser::getCSF(std::istream& is)
}
getwords(currentWords,is);
} else {
if(NTOT <= 50)
if(NTOT < 50)
getwords(currentWords,is);
else {
getwords(currentWords,is);
@ -1102,7 +1126,7 @@ void GamesAsciiParser::getORMAS(std::istream& is)
abort();
}
parsewords(aline.c_str(),currentWords);
ci_nea = atoi(currentWords[4].c_str());
ci_nca = ci_ncb = atoi(currentWords[4].c_str());
if(!lookFor(is,"NUMBER OF ACTIVE ORBITALS",aline))
{
@ -1119,7 +1143,8 @@ void GamesAsciiParser::getORMAS(std::istream& is)
abort();
}
parsewords(aline.c_str(),currentWords);
ci_nea = atoi(currentWords[4].c_str());
//ci_nea = atoi(currentWords[4].c_str());
ci_nea = atoi(currentWords[6].c_str());
if(!lookFor(is,"NUMBER OF BETA ELECTRONS",aline))
{
@ -1127,7 +1152,14 @@ void GamesAsciiParser::getORMAS(std::istream& is)
abort();
}
parsewords(aline.c_str(),currentWords);
ci_neb = atoi(currentWords[4].c_str());
//ci_neb = atoi(currentWords[4].c_str());
ci_neb = atoi(currentWords[6].c_str());
cout <<"ORMAS: nea,neb,ncore,nact: "
<<ci_nea <<" "
<<ci_neb <<" "
<<ci_nca <<" "
<<nactive <<"\n";
int ds=SpinMultiplicity-1;
int neb= (NumberOfEls-ds)/2;

4
src/QMCTools/QMCGaussianParserBase.cpp
View File

@ -578,12 +578,13 @@ QMCGaussianParserBase::createMultiDeterminantSet() {
ci_size=0;
for(int i=0; i<CIcoeff.size(); i++) if(fabs(CIcoeff[i]) > ci_threshold) ci_size++;
std::ostringstream nstates,cisize,cinca,cincb,cinea,cineb;
std::ostringstream nstates,cisize,cinca,cincb,cinea,cineb,ci_thr;
cisize <<ci_size; nstates <<ci_nstates;
cinca <<ci_nca;
cincb <<ci_ncb;
cinea <<ci_nea;
cineb <<ci_neb;
ci_thr <<ci_threshold;
xmlNewProp(detlist,(const xmlChar*)"size",(const xmlChar*)cisize.str().c_str());
xmlNewProp(detlist,(const xmlChar*)"type",(const xmlChar*)"DETS");
@ -592,6 +593,7 @@ QMCGaussianParserBase::createMultiDeterminantSet() {
xmlNewProp(detlist,(const xmlChar*)"nea",(const xmlChar*)cinea.str().c_str());
xmlNewProp(detlist,(const xmlChar*)"neb",(const xmlChar*)cineb.str().c_str());
xmlNewProp(detlist,(const xmlChar*)"nstates",(const xmlChar*)nstates.str().c_str());
xmlNewProp(detlist,(const xmlChar*)"cutoff",(const xmlChar*)ci_thr.str().c_str());
if(CIcoeff.size() == 0) {
cerr<<" CI configuration list is empty. \n";

32
src/QMCWaveFunctions/Experimental/LCOrbitalSetWithCorrection.cpp
View File

@ -43,9 +43,11 @@ namespace qmcplusplus {
{
int bss = myBasisSet->LOBasis[i]->BasisSetSize;
for(int k=0; k<bss; k++)
rmv[cnt++] = (myBasisSet->LOBasis[i]->RnlID[myBasisSet->LOBasis[i]->NL[k]])[1]==0;
rmv[cnt++] = (((myBasisSet->LOBasis[i]->RnlID[myBasisSet->LOBasis[i]->NL[k]])[1]==0)&&(corrCenter[i]));
}
BS* Eta;
/* not used any more, restore if you want to test and propagate changes due to corrCenter
BS* Eta = new BS(*myBasisSet);
for(int i=0; i<nUniqCenters; i++)
{
@ -119,7 +121,7 @@ namespace qmcplusplus {
// now set basis set sizes
Eta->setBasisSetSize(-10);
Eta->resize(Eta->NumTargets);
*/
return Eta;
}
@ -141,7 +143,7 @@ namespace qmcplusplus {
// WARNING, assuming that COT has to be NGOrbital, otherwise problems
COT* cur = (COT*) myBasisSet->LOBasis[i];
for(int k=0; k<bss; k++)
rmv[cnt++] = (cur->RnlID[cur->NL[k]])[1]==0;
rmv[cnt++] = (((cur->RnlID[cur->NL[k]])[1]==0)&&(corrCenter[i]));;
}
else
{
@ -297,7 +299,17 @@ app_log() <<"cuspFile: " <<cuspInfoFile <<endl;
for(int iat=0; iat<Z.size();iat++) {
Z[iat] = tspecies(iz,sourcePtcl->GroupID[iat]);
}
int numAttrib = tspecies.numAttributes();
int cctag = tspecies.addAttribute("cuspCorr");
corrCenter.resize(Z.size(),"true");
if(cctag < numAttrib) { // parameter is not new
for(int iat=0; iat<Z.size();iat++) {
if( tspecies(cctag,sourcePtcl->GroupID[iat]) != 0) {
corrCenter[iat] = false;
app_log() <<"Not using cusp correction algorithm in atoms of type: " <<tspecies.speciesName[sourcePtcl->GroupID[iat]] <<endl;
}
}
}
if(Rcut < 0) Rcut=0.1;
//int indexRc=-1;
@ -381,10 +393,12 @@ app_log() <<"cuspFile: " <<cuspInfoFile <<endl;
num0<<k;
xmlNewProp(orb,(const xmlChar*)"num",(const xmlChar*)num0.str().c_str());
bool corrO = false;
for(int ip=0; ip<dummyLO1->C.cols(); ip++) {
if(std::fabs(dummyLO1->C(k,ip)) > 1e-8) {
corrO = true;
break;
if(corrCenter[i]) {
for(int ip=0; ip<dummyLO1->C.cols(); ip++) {
if(std::fabs(dummyLO1->C(k,ip)) > 1e-8) {
corrO = true;
break;
}
}
}
if(corrO) {
@ -475,7 +489,7 @@ app_log() <<"cuspFile: " <<cuspInfoFile <<endl;
corrBasisSet->setBasisSetSize(-1);
BS *dum3 = extractHighYLM(rmv);
dum3->setBasisSetSize(-1);
//dum3->setBasisSetSize(-1);
int norb=0, cnt=0;
for(int i=0; i<rmv.size(); i++) if(!rmv[i]) norb++;
app_log()<<"Found " <<rmv.size()-norb <<" spherically symmetric basis functions and " <<norb <<" non symmetric ones. \n";

1
src/QMCWaveFunctions/Experimental/LCOrbitalSetWithCorrection.h
View File

@ -203,6 +203,7 @@ namespace qmcplusplus {
///pointer to the basis set
BS* myBasisSet;
CorrectingOrbitalBasisSet<COT>* corrBasisSet;
vector<bool> corrCenter;
///target ParticleSet
ParticleSet* targetPtcl;

53
src/QMCWaveFunctions/Fermion/MultiDiracDeterminantBase.cpp
View File

@ -44,6 +44,8 @@ namespace qmcplusplus {
FirstIndex = first;
DetCalculator.resize(nel);
resize(nel,norb);
// mmorales; remove later
// testDets();
}
void MultiDiracDeterminantBase::createDetData(ci_configuration2& ref, vector<int>& data,
@ -95,13 +97,25 @@ namespace qmcplusplus {
//{ cout<<"MDD: " <<str <<" " <<n <<endl; cout.flush(); }
void MultiDiracDeterminantBase::evaluateForWalkerMove(ParticleSet& P)
void MultiDiracDeterminantBase::evaluateForWalkerMove(ParticleSet& P, bool fromScratch)
{
evalWTimer.start();
Phi->evaluate_notranspose(P,FirstIndex,LastIndex,psiM,dpsiM,d2psiM);
if(fromScratch) Phi->evaluate_notranspose(P,FirstIndex,LastIndex,psiM,dpsiM,d2psiM);
if(NumPtcls==1) {
APP_ABORT("Evaluate Log with 1 particle in MultiDiracDeterminantBase is potentially dangerous. Fix later");
//APP_ABORT("Evaluate Log with 1 particle in MultiDiracDeterminantBase is potentially dangerous. Fix later");
vector<ci_configuration2>::iterator it(confgList.begin());
vector<ci_configuration2>::iterator last(confgList.end());
ValueVector_t::iterator det(detValues.begin());
ValueMatrix_t::iterator lap(lapls.begin());
GradMatrix_t::iterator grad(grads.begin());
while(it != last) {
int orb = (it++)->occup[0];
*(det++) = psiM(0,orb);
*(lap++) = d2psiM(0,orb);
*(grad++) = dpsiM(0,orb);
}
} else {
InverseTimer.start();
@ -164,11 +178,6 @@ namespace qmcplusplus {
} // NumPtcls==1
psiMinv_temp = psiMinv;
// only used with ORB_PBYP_ALL,
// do I need this here???
//psiM_temp = psiM;
//dpsiM_temp = dpsiM;
//d2psiM_temp = d2psiM;
evalWTimer.stop();
}
@ -178,8 +187,7 @@ namespace qmcplusplus {
PooledData<RealType>& buf, bool fromscratch)
{
if(fromscratch)
evaluateForWalkerMove(P);
evaluateForWalkerMove(P,(fromscratch || UpdateMode == ORB_PBYP_RATIO) );
buf.put(psiM.first_address(),psiM.last_address());
buf.put(FirstAddressOfdpsiM,LastAddressOfdpsiM);
@ -225,6 +233,7 @@ namespace qmcplusplus {
*/
void MultiDiracDeterminantBase::acceptMove(ParticleSet& P, int iat)
{
WorkingIndex = iat-FirstIndex;
switch(UpdateMode)
{
@ -245,6 +254,7 @@ namespace qmcplusplus {
// grads(i,WorkingIndex) = new_grads(i,WorkingIndex);
std::copy(psiV.begin(),psiV.end(),psiM[WorkingIndex]);
std::copy(dpsiV.begin(),dpsiV.end(),dpsiM[WorkingIndex]);
std::copy(d2psiV.begin(),d2psiV.end(),d2psiM[WorkingIndex]);
break;
default:
psiMinv = psiMinv_temp;
@ -264,7 +274,6 @@ namespace qmcplusplus {
/** move was rejected. copy the real container to the temporary to move on
*/
void MultiDiracDeterminantBase::restore(int iat) {
WorkingIndex = iat-FirstIndex;
psiMinv_temp = psiMinv;
for(int i=0; i<NumOrbitals; i++)
@ -320,7 +329,11 @@ namespace qmcplusplus {
buildTableTimer("MultiDiracDeterminantBase::buildTable"),
evalWTimer("MultiDiracDeterminantBase::evalW"),
evalOrbTimer("MultiDiracDeterminantBase::evalOrb"),
readMatTimer("MultiDiracDeterminantBase::readMat")
evalOrb1Timer("MultiDiracDeterminantBase::evalOrbGrad"),
readMatTimer("MultiDiracDeterminantBase::readMat"),
readMatGradTimer("MultiDiracDeterminantBase::readMatGrad"),
buildTableGradTimer("MultiDiracDeterminantBase::buildTableGrad"),
ExtraStuffTimer("MultiDiracDeterminantBase::ExtraStuff")
{
setDetInfo(s.ReferenceDeterminant,s.confgList);
registerTimers();
@ -352,7 +365,11 @@ namespace qmcplusplus {
buildTableTimer("MultiDiracDeterminantBase::buildTable"),
evalWTimer("MultiDiracDeterminantBase::evalW"),
evalOrbTimer("MultiDiracDeterminantBase::evalOrb"),
readMatTimer("MultiDiracDeterminantBase::readMat")
evalOrb1Timer("MultiDiracDeterminantBase::evalOrbGrad"),
readMatTimer("MultiDiracDeterminantBase::readMat"),
readMatGradTimer("MultiDiracDeterminantBase::readMatGrad"),
buildTableGradTimer("MultiDiracDeterminantBase::buildTableGrad"),
ExtraStuffTimer("MultiDiracDeterminantBase::ExtraStuff")
{
Optimizable=true;
OrbitalName="MultiDiracDeterminantBase";
@ -385,7 +402,7 @@ namespace qmcplusplus {
LastAddressOfdpsiM = FirstAddressOfdpsiM + NumPtcls*NumOrbitals*DIM;
}
evaluateForWalkerMove(P);
evaluateForWalkerMove(P,true);
//add the data:
buf.add(psiM.first_address(),psiM.last_address());
@ -492,7 +509,11 @@ namespace qmcplusplus {
buildTableTimer.reset();
readMatTimer.reset();
evalOrbTimer.reset();
evalOrb1Timer.reset();
evalWTimer.reset();
ExtraStuffTimer.reset();
buildTableGradTimer.reset();
readMatGradTimer.reset();
TimerManager.addTimer (&UpdateTimer);
TimerManager.addTimer (&RatioTimer);
TimerManager.addTimer (&InverseTimer);
@ -500,6 +521,10 @@ namespace qmcplusplus {
TimerManager.addTimer (&readMatTimer);
TimerManager.addTimer (&evalWTimer);
TimerManager.addTimer (&evalOrbTimer);
TimerManager.addTimer (&evalOrb1Timer);
TimerManager.addTimer (&buildTableGradTimer);
TimerManager.addTimer (&readMatGradTimer);
TimerManager.addTimer (&ExtraStuffTimer);
}

488
src/QMCWaveFunctions/Fermion/MultiDiracDeterminantBase.h
View File

@ -34,7 +34,8 @@ namespace qmcplusplus {
public:
bool Optimizable;
void registerTimers();
NewTimer UpdateTimer, RatioTimer, InverseTimer, buildTableTimer, readMatTimer, evalWTimer, evalOrbTimer;
NewTimer UpdateTimer, RatioTimer, InverseTimer, buildTableTimer, readMatTimer, evalWTimer, evalOrbTimer,evalOrb1Timer;
NewTimer readMatGradTimer,buildTableGradTimer,ExtraStuffTimer;
// Optimizable parameters
opt_variables_type myVars;
@ -215,6 +216,101 @@ namespace qmcplusplus {
// this works with confgList, which shouldn't change during a simulation
void createDetData(ci_configuration2& ref, vector<int>& data,
vector<pair<int,int> >& pairs, vector<double>& sign);
/*
void testDets() {
vector<int> indx(16);
ValueMatrix_t Mat(8,8);
NewTimer dummy("dummy");
dummy.reset();
ifstream in("matrix.txt");
for(int i=0; i<8; i++)
for(int j=0; j<8; j++)
in>>Mat(i,j);
in.close();
vector<double> dets(9);
dets[2] = -0.085167; // 2x2
dets[3] = 0.051971; // 3x3
dets[4] = 0.029865; // 4x4
dets[5] = -0.055246; // 5x5
dets[6] = -0.0087414; // 6x6
dets[7] = 0.025257; // 7x7
dets[8] = 0.067119; // 8x8
for(int q=2; q<=8; q++) {
int k=0;
for(int i=0; i<q; i++)
indx[k++] = i;
for(int i=0; i<q; i++)
indx[k++] = i;
cout<<"determinant of " <<q <<": " <<dets[q] <<" - " <<CalculateRatioFromMatrixElements(q,Mat,indx.begin()) <<endl;
}
int k=0;
for(int i=0; i<4; i++)
indx[k++] = i;
for(int i=0; i<4; i++)
indx[k++] = i;
cout<<"Timing n=4 \n";
ValueType res;
dummy.start();
for(int i=0; i<1000; i++)
res=CalculateRatioFromMatrixElements(4,Mat,indx.begin());
dummy.stop();
double direct = dummy.get_total();
dummy.reset();
dummy.start();
for(int i=0; i<1000; i++)
res=NewCalculateRatioFromMatrixElements(4,Mat,indx.begin());
dummy.stop();
double func = dummy.get_total();
cout<<"Direct, function: " <<direct <<" " <<func <<endl <<endl;
k=0;
for(int i=0; i<5; i++)
indx[k++] = i;
for(int i=0; i<5; i++)
indx[k++] = i;
cout<<"Timing n=5 \n";
dummy.start();
for(int i=0; i<1000; i++)
res=CalculateRatioFromMatrixElements(5,Mat,indx.begin());
dummy.stop();
direct = dummy.get_total();
dummy.reset();
dummy.start();
for(int i=0; i<1000; i++)
res=NewCalculateRatioFromMatrixElements(5,Mat,indx.begin());
dummy.stop();
func = dummy.get_total();
cout<<"Direct, function: " <<direct <<" " <<func <<endl <<endl;
k=0;
for(int i=0; i<6; i++)
indx[k++] = i;
for(int i=0; i<6; i++)
indx[k++] = i;
cout<<"Timing n=6 \n";
dummy.start();
for(int i=0; i<1000; i++)
res=CalculateRatioFromMatrixElements(6,Mat,indx.begin());
dummy.stop();
direct = dummy.get_total();
dummy.reset();
dummy.start();
for(int i=0; i<1000; i++)
res=NewCalculateRatioFromMatrixElements(6,Mat,indx.begin());
dummy.stop();
func = dummy.get_total();
cout<<"Direct, function: " <<direct <<" " <<func <<endl <<endl;
exit(1);
}
*/
inline ValueType CalculateRatioFromMatrixElements(int n, ValueMatrix_t& dotProducts, vector<int>::iterator it)
{
@ -285,6 +381,7 @@ namespace qmcplusplus {
dotProducts(i5,a1),dotProducts(i5,a2),dotProducts(i5,a3),dotProducts(i5,a4),dotProducts(i5,a5));
break;
}
/*
case 6:
{
register int i1=*it;
@ -308,12 +405,39 @@ namespace qmcplusplus {
dotProducts(i6,a1),dotProducts(i6,a2),dotProducts(i6,a3),dotProducts(i6,a4),dotProducts(i6,a5),dotProducts(i6,a6));
break;
}
*/
default:
{
return DetCalculator.evaluate(dotProducts,it,n);
}
}
}
}
/*
inline ValueType NewCalculateRatioFromMatrixElements(int n, ValueMatrix_t& dotProducts, vector<int>::iterator it)
{
switch(n)
{
case 0:
return 1.0;
case 1:
return dotProducts(*it,*(it+1));
break;
case 2:
{
register int i=*it;
register int j=*(it+1);
register int a=*(it+2);
register int b=*(it+3);
return dotProducts(i,a)*dotProducts(j,b)-dotProducts(i,b)*dotProducts(j,a);
break;
}
default:
{
return DetCalculator.evaluate(dotProducts,it,n);
}
}
}
*/
inline void BuildDotProductsAndCalculateRatios(int ref, int iat, ValueVector_t& ratios, ValueMatrix_t &psiinv, ValueMatrix_t &psi, ValueMatrix_t& dotProducts, vector<int>& data, vector<pair<int,int> >& pairs, vector<double>& sign)
{
@ -354,13 +478,16 @@ namespace qmcplusplus {
{
ValueType det0 = ratios(ref,iat)[dx];
buildTableGradTimer.start();
int num=psi.extent(1);
vector<pair<int,int> >::iterator it(pairs.begin()), last(pairs.end());
while(it != last) {
dotProducts((*it).first,(*it).second) = dot(psiinv[(*it).first],psi[(*it).second],num);
it++;
}
buildTableGradTimer.stop();
readMatGradTimer.start();
vector<int>::iterator it2 = data.begin();
int count= 0; // number of determinants processed
while(it2 != data.end()) {
@ -374,6 +501,7 @@ namespace qmcplusplus {
count++;
it2+=3*n+1;
}
readMatGradTimer.stop();
}
inline void BuildDotProductsAndCalculateRatios(int ref, int iat, ValueMatrix_t& ratios, ValueMatrix_t& psiinv, ValueMatrix_t& psi, ValueMatrix_t& dotProducts, vector<int>& data, vector<pair<int,int> >& pairs, vector<double>& sign)
@ -440,68 +568,100 @@ namespace qmcplusplus {
evalOrbTimer.stop();
WorkingIndex = iat-FirstIndex;
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
if(NumPtcls==1) {
vector<ci_configuration2>::iterator it(confgList.begin());
vector<ci_configuration2>::iterator last(confgList.end());
ValueVector_t::iterator det(new_detValues.begin());
while(it != last) {
int orb = (it++)->occup[0];
*(det++) = psiV(orb);
}
} else {
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
// mmorales: the only reason this is here is because
// NonlocalECP do not necessarily call rejectMove after
// calling ratio(), and even if the move is rejected
// this matrix needs to be restored
// If we always restore after ratio, then this is not needed
// For efficiency reasons, I don't do this for ratioGrad or ratio(P,dG,dL)
psiMinv_temp = psiMinv;
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = psiV(*(it++));
ValueType ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
new_detValues[ReferenceDeterminant] = ratioRef * detValues[ReferenceDeterminant];
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiV(i);
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
ExtraStuffTimer.start();
psiMinv_temp = psiMinv;
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = psiV(*(it++));
ValueType ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
new_detValues[ReferenceDeterminant] = ratioRef * detValues[ReferenceDeterminant];
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiV(i);
ExtraStuffTimer.stop();
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
// check comment above
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
}
RatioTimer.stop();
}
inline void
evaluateDetsAndGradsForPtclMove(ParticleSet& P, int iat) {
UpdateMode=ORB_PBYP_PARTIAL;
evalOrb1Timer.start();
Phi->evaluate(P,iat,psiV,dpsiV,d2psiV);
//mmorales: check comment above
psiMinv_temp = psiMinv;
evalOrb1Timer.stop();
WorkingIndex = iat-FirstIndex;
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
GradType ratioGradRef;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = psiV(*it);
ratioGradRef += psiMinv_temp(i,WorkingIndex)*dpsiV(*it);
it++;
}
ValueType ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
new_grads(ReferenceDeterminant,WorkingIndex) = ratioGradRef*detValues[ReferenceDeterminant];
new_detValues[ReferenceDeterminant] = ratioRef * detValues[ReferenceDeterminant];
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiV[i];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
if(NumPtcls==1) {
for(int idim=0; idim<3; idim++) {
//dpsiMinv = psiMinv_temp;
dpsiMinv = psiMinv;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiV(*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,WorkingIndex,ratioGradRef[idim]);
vector<ci_configuration2>::iterator it(confgList.begin());
vector<ci_configuration2>::iterator last(confgList.end());
ValueVector_t::iterator det(new_detValues.begin());
GradMatrix_t::iterator grad(new_grads.begin());
while(it != last) {
int orb = (it++)->occup[0];
*(det++) = psiV[orb];
*(grad++) = dpsiV[orb];
}
} else {
ExtraStuffTimer.start();
//mmorales: check comment above
psiMinv_temp = psiMinv;
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
GradType ratioGradRef;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = psiV(*it);
ratioGradRef += psiMinv_temp(i,WorkingIndex)*dpsiV(*it);
it++;
}
ValueType ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
new_grads(ReferenceDeterminant,WorkingIndex) = ratioGradRef*detValues[ReferenceDeterminant];
new_detValues[ReferenceDeterminant] = ratioRef * detValues[ReferenceDeterminant];
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = dpsiV[i][idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
TpsiM(i,WorkingIndex) = psiV[i];
ExtraStuffTimer.stop();
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int idim=0; idim<3; idim++) {
ExtraStuffTimer.start();
//dpsiMinv = psiMinv_temp;
dpsiMinv = psiMinv;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiV(*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,WorkingIndex,ratioGradRef[idim]);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = dpsiV[i][idim];
ExtraStuffTimer.stop();
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
// check comment above
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
}
}
inline void
@ -510,27 +670,36 @@ namespace qmcplusplus {
WorkingIndex = iat-FirstIndex;
vector<int>::iterator it;
for(int idim=0; idim<3; idim++) {
//dpsiMinv = psiMinv_temp;
dpsiMinv = psiMinv;
it = confgList[ReferenceDeterminant].occup.begin();
ValueType ratioG = 0.0;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = dpsiM(WorkingIndex,*it)[idim];
ratioG += psiMinv(i,WorkingIndex)*dpsiM(WorkingIndex,*it)[idim];
it++;
if(NumPtcls==1) {
vector<ci_configuration2>::iterator it(confgList.begin());
vector<ci_configuration2>::iterator last(confgList.end());
GradMatrix_t::iterator grad(grads.begin());
while(it != last) {
int orb = (it++)->occup[0];
*(grad++) = dpsiM(0,orb);
}
} else {
for(int idim=0; idim<3; idim++) {
//dpsiMinv = psiMinv_temp;
dpsiMinv = psiMinv;
it = confgList[ReferenceDeterminant].occup.begin();
ValueType ratioG = 0.0;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = dpsiM(WorkingIndex,*it)[idim];
ratioG += psiMinv(i,WorkingIndex)*dpsiM(WorkingIndex,*it)[idim];
it++;
}
grads(ReferenceDeterminant,WorkingIndex)[idim] = ratioG*detValues[ReferenceDeterminant];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,WorkingIndex,ratioG);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = dpsiM(WorkingIndex,i)[idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
grads(ReferenceDeterminant,WorkingIndex)[idim] = ratioG*detValues[ReferenceDeterminant];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,WorkingIndex,ratioG);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = dpsiM(WorkingIndex,i)[idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
// check comment above
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
}
}
inline void
@ -538,105 +707,120 @@ namespace qmcplusplus {
UpdateMode=ORB_PBYP_ALL;
Phi->evaluate(P,iat,psiV,dpsiV,d2psiV);
//mmorales: check comment above
psiMinv_temp = psiMinv;
WorkingIndex = iat-FirstIndex;
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
//GradType ratioGradRef;
//ValueType ratioLapl = 0.0;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = psiV(*it);
//ratioGradRef += psiMinv_temp(i,WorkingIndex)*dpsiV(*it);
//ratioLapl += psiMinv_temp(i,WorkingIndex)*d2psiV(*it);
it++;
}
ValueType det0, ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
//new_lapls(ReferenceDeterminant,WorkingIndex) = ratioLapl*detValues[ReferenceDeterminant];
//new_grads(ReferenceDeterminant,WorkingIndex) = ratioGradRef*detValues[ReferenceDeterminant];
det0 = ratioRef * detValues[ReferenceDeterminant];
new_detValues[ReferenceDeterminant] = det0;
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiV[i];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int jat=0; jat<NumPtcls; jat++) {
it = confgList[ReferenceDeterminant].occup.begin();
GradType gradRatio = 0.0;
ValueType ratioLapl = 0.0;
if(jat==WorkingIndex) {
for(int i=0; i<NumPtcls; i++) {
gradRatio += psiMinv_temp(i,jat)*dpsiV(*it);
ratioLapl += psiMinv_temp(i,jat)*d2psiV(*it);
it++;
}
new_grads(ReferenceDeterminant,jat) = det0*gradRatio;
new_lapls(ReferenceDeterminant,jat) = det0*ratioLapl;
for(int idim=0; idim<3; idim++) {
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiV(*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,gradRatio[idim]);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = dpsiV[i][idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = d2psiV(*(it++));
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,ratioLapl);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = d2psiV[i];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_lapls,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = psiV[i];
} else {
for(int i=0; i<NumPtcls; i++) {
gradRatio += psiMinv_temp(i,jat)*dpsiM(jat,*it);
ratioLapl += psiMinv_temp(i,jat)*d2psiM(jat,*it);
it++;
}
new_grads(ReferenceDeterminant,jat) = det0*gradRatio;
new_lapls(ReferenceDeterminant,jat) = det0*ratioLapl;
for(int idim=0; idim<3; idim++) {
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiM(jat,*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,gradRatio[idim]);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = dpsiM(jat,i)[idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = d2psiM(jat,*(it++));
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,ratioLapl);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = d2psiM(jat,i);
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_lapls,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = psiM(jat,i);
if(NumPtcls==1) {
vector<ci_configuration2>::iterator it(confgList.begin());
vector<ci_configuration2>::iterator last(confgList.end());
ValueVector_t::iterator det(new_detValues.begin());
ValueMatrix_t::iterator lap(new_lapls.begin());
GradMatrix_t::iterator grad(new_grads.begin());
while(it != last) {
int orb = (it++)->occup[0];
*(det++) = psiV(orb);
*(lap++) = d2psiV(orb);
*(grad++) = dpsiV(orb);
}
} // jat
} else {
//mmorales: check comment above
psiMinv_temp = psiMinv;
vector<int>::iterator it(confgList[ReferenceDeterminant].occup.begin());
//GradType ratioGradRef;
//ValueType ratioLapl = 0.0;
for(int i=0; i<NumPtcls; i++) {
psiV_temp[i] = psiV(*it);
//ratioGradRef += psiMinv_temp(i,WorkingIndex)*dpsiV(*it);
//ratioLapl += psiMinv_temp(i,WorkingIndex)*d2psiV(*it);
it++;
}
ValueType det0, ratioRef = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
//new_lapls(ReferenceDeterminant,WorkingIndex) = ratioLapl*detValues[ReferenceDeterminant];
//new_grads(ReferenceDeterminant,WorkingIndex) = ratioGradRef*detValues[ReferenceDeterminant];
det0 = ratioRef * detValues[ReferenceDeterminant];
new_detValues[ReferenceDeterminant] = det0;
InverseUpdateByColumn(psiMinv_temp,psiV_temp,workV1,workV2,WorkingIndex,ratioRef);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiV[i];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,WorkingIndex,new_detValues,psiMinv_temp,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int jat=0; jat<NumPtcls; jat++) {
it = confgList[ReferenceDeterminant].occup.begin();
GradType gradRatio = 0.0;
ValueType ratioLapl = 0.0;
if(jat==WorkingIndex) {
for(int i=0; i<NumPtcls; i++) {
gradRatio += psiMinv_temp(i,jat)*dpsiV(*it);
ratioLapl += psiMinv_temp(i,jat)*d2psiV(*it);
it++;
}
new_grads(ReferenceDeterminant,jat) = det0*gradRatio;
new_lapls(ReferenceDeterminant,jat) = det0*ratioLapl;
for(int idim=0; idim<3; idim++) {
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiV(*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,gradRatio[idim]);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = dpsiV[i][idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = d2psiV(*(it++));
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,ratioLapl);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = d2psiV[i];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_lapls,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = psiV[i];
} else {
for(int i=0; i<NumPtcls; i++) {
gradRatio += psiMinv_temp(i,jat)*dpsiM(jat,*it);
ratioLapl += psiMinv_temp(i,jat)*d2psiM(jat,*it);
it++;
}
new_grads(ReferenceDeterminant,jat) = det0*gradRatio;
new_lapls(ReferenceDeterminant,jat) = det0*ratioLapl;
for(int idim=0; idim<3; idim++) {
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = dpsiM(jat,*(it++))[idim];
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,gradRatio[idim]);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = dpsiM(jat,i)[idim];
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_grads,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns,idim);
}
dpsiMinv = psiMinv_temp;
it = confgList[ReferenceDeterminant].occup.begin();
for(int i=0; i<NumPtcls; i++)
psiV_temp[i] = d2psiM(jat,*(it++));
InverseUpdateByColumn(dpsiMinv,psiV_temp,workV1,workV2,jat,ratioLapl);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = d2psiM(jat,i);
BuildDotProductsAndCalculateRatios(ReferenceDeterminant,jat,new_lapls,dpsiMinv,TpsiM,dotProducts,detData,uniquePairs,DetSigns);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,jat) = psiM(jat,i);
}
} // jat
// check comment above
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
for(int i=0; i<NumOrbitals; i++)
TpsiM(i,WorkingIndex) = psiM(WorkingIndex,i);
}
}
// full evaluation of all the structures from scratch, used in evaluateLog for example
void evaluateForWalkerMove(ParticleSet& P);
void evaluateForWalkerMove(ParticleSet& P, bool fromScratch=true);
///total number of particles
int NP;

9
src/QMCWaveFunctions/Fermion/MultiDiracDeterminantBase_help.h
View File

@ -66,13 +66,14 @@ namespace qmcplusplus {
T a51, T a52, T a53, T a54, T a55)
{
//return a11*evaluate(a22,a23,a24,a25,a32,a33,a34,a35,a42,a43,a44,a45,a52,a53,a54,a55)
- a21*evaluate(a12,a13,a14,a15,a32,a33,a34,a35,a42,a43,a44,a45,a52,a53,a54,a55)
+ a31*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a42,a43,a44,a45,a52,a53,a54,a55)
- a41*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a32,a33,a34,a35,a52,a53,a54,a55)
+ a51*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a32,a33,a34,a35,a42,a43,a44,a45);
// - a21*evaluate(a12,a13,a14,a15,a32,a33,a34,a35,a42,a43,a44,a45,a52,a53,a54,a55)
// + a31*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a42,a43,a44,a45,a52,a53,a54,a55)
// - a41*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a32,a33,a34,a35,a52,a53,a54,a55)
// + a51*evaluate(a12,a13,a14,a15,a22,a23,a24,a25,a32,a33,a34,a35,a42,a43,a44,a45);
return (a11*(a22*(a33*(a44*a55-a54*a45)-a43*(a34*a55-a54*a35)+a53*(a34*a45-a44*a35))-a32*(a23*(a44*a55-a54*a45)-a43*(a24*a55-a54*a25)+a53*(a24*a45-a44*a25))+a42*(a23*(a34*a55-a54*a35)-a33*(a24*a55-a54*a25)+a53*(a24*a35-a34*a25))-a52*(a23*(a34*a45-a44*a35)-a33*(a24*a45-a44*a25)+a43*(a24*a35-a34*a25)))-a21*(a12*(a33*(a44*a55-a54*a45)-a43*(a34*a55-a54*a35)+a53*(a34*a45-a44*a35))-a32*(a13*(a44*a55-a54*a45)-a43*(a14*a55-a54*a15)+a53*(a14*a45-a44*a15))+a42*(a13*(a34*a55-a54*a35)-a33*(a14*a55-a54*a15)+a53*(a14*a35-a34*a15))-a52*(a13*(a34*a45-a44*a35)-a33*(a14*a45-a44*a15)+a43*(a14*a35-a34*a15)))+a31*(a12*(a23*(a44*a55-a54*a45)-a43*(a24*a55-a54*a25)+a53*(a24*a45-a44*a25))-a22*(a13*(a44*a55-a54*a45)-a43*(a14*a55-a54*a15)+a53*(a14*a45-a44*a15))+a42*(a13*(a24*a55-a54*a25)-a23*(a14*a55-a54*a15)+a53*(a14*a25-a24*a15))-a52*(a13*(a24*a45-a44*a25)-a23*(a14*a45-a44*a15)+a43*(a14*a25-a24*a15)))-a41*(a12*(a23*(a34*a55-a54*a35)-a33*(a24*a55-a54*a25)+a53*(a24*a35-a34*a25))-a22*(a13*(a34*a55-a54*a35)-a33*(a14*a55-a54*a15)+a53*(a14*a35-a34*a15))+a32*(a13*(a24*a55-a54*a25)-a23*(a14*a55-a54*a15)+a53*(a14*a25-a24*a15))-a52*(a13*(a24*a35-a34*a25)-a23*(a14*a35-a34*a15)+a33*(a14*a25-a24*a15)))+a51*(a12*(a23*(a34*a45-a44*a35)-a33*(a24*a45-a44*a25)+a43*(a24*a35-a34*a25))-a22*(a13*(a34*a45-a44*a35)-a33*(a14*a45-a44*a15)+a43*(a14*a35-a34*a15))+a32*(a13*(a24*a45-a44*a25)-a23*(a14*a45-a44*a15)+a43*(a14*a25-a24*a15))-a42*(a13*(a24*a35-a34*a25)-a23*(a14*a35-a34*a15)+a33*(a14*a25-a24*a15))));
}
// mmorales: THERE IS SOMETHING WRONG WITH THIS ROUTINE, BUT I DON'T USE IT ANYMORE !!!
inline T evaluate(T a11, T a12, T a13, T a14, T a15, T a16,
T a21, T a22, T a23, T a24, T a25, T a26,
T a31, T a32, T a33, T a34, T a35, T a36,

209
src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp
View File

@ -26,12 +26,15 @@ namespace qmcplusplus {
RatioAllTimer("MultiSlaterDeterminant::ratio(all)"),
Ratio1AllTimer("MultiSlaterDeterminant::detEval_ratio(all)"),
UpdateTimer("MultiSlaterDeterminant::updateBuffer"),
EvaluateTimer("MultiSlaterDeterminant::evaluate")
AccRejTimer("MultiSlaterDeterminant::Accept_Reject"),
EvaluateTimer("MultiSlaterDeterminant::evaluate"),
evalOrbTimer("MultiSlaterDeterminant::evalOrbGrad")
{
registerTimers();
//Optimizable=true;
Optimizable=false;
OrbitalName="MultiSlaterDeterminant";
usingCSF=false;
FirstIndex_up = targetPtcl.first(0);
LastIndex_up = targetPtcl.last(0);
@ -59,6 +62,12 @@ namespace qmcplusplus {
clone->C2node_up=C2node_up;
clone->C2node_dn=C2node_dn;
clone->resize(dets_up.size(),dets_dn.size());
if (usingCSF)
{
clone->CSFcoeff=CSFcoeff;
clone->CSFexpansion=CSFexpansion;
clone->DetsPerCSF=DetsPerCSF;
}
// SPOSetProxyForMSD* spo = clone->spo_up;
// spo->occup.resize(uniqueConfg_up.size(),clone->nels_up);
@ -252,7 +261,9 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
UpdateMode=ORB_PBYP_PARTIAL;
if(DetID[iat] == 0) {
RatioGradTimer.start();
evalOrbTimer.start();
spo_up->evaluateAllForPtclMove(P,iat);
evalOrbTimer.stop();
Ratio1GradTimer.start();
grad_temp=0.0;
@ -279,7 +290,9 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
return curRatio;
} else {
RatioGradTimer.start();
evalOrbTimer.start();
spo_dn->evaluateAllForPtclMove(P,iat);
evalOrbTimer.stop();
Ratio1GradTimer.start();
grad_temp=0.0;
@ -437,48 +450,47 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
{
UpdateMode=ORB_PBYP_RATIO;
if(DetID[iat] == 0) {
RatioTimer.start();
RatioTimer.start();
spo_up->evaluateForPtclMove(P,iat);
Ratio1Timer.start();
Ratio1Timer.start();
for(int i=0; i<dets_up.size(); i++) {
spo_up->prepareFor(i);
detsRatios[i]=dets_up[i]->ratio(P,iat);
}
Ratio1Timer.stop();
Ratio1Timer.stop();
vector<int>::iterator upC(C2node_up.begin()),dnC(C2node_dn.begin());
vector<RealType>::iterator it(C.begin()),last(C.end());
ValueType psiOld=0.0,psiNew=0.0;
for(int i=0; i<C.size(); i++){
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType tmp2 = C[i]*detValues_up[upC]*detValues_dn[dnC];
psiOld += tmp2;
psiNew += tmp2*detsRatios[upC];
while(it != last) {
ValueType tmp = (*it)*detValues_up[*upC]*detValues_dn[*dnC];
psiNew += tmp*detsRatios[*upC];
psiOld += tmp;
it++;upC++;dnC++;
}
curRatio = psiNew/psiOld;
RatioTimer.stop();
RatioTimer.stop();
return curRatio;
} else {
RatioTimer.start();
RatioTimer.start();
spo_dn->evaluateForPtclMove(P,iat);
Ratio1Timer.start();
Ratio1Timer.start();
for(int i=0; i<dets_dn.size(); i++) {
spo_dn->prepareFor(i);
detsRatios[i]=dets_dn[i]->ratio(P,iat);
}
Ratio1Timer.stop();
Ratio1Timer.stop();
vector<int>::iterator upC(C2node_up.begin()),dnC(C2node_dn.begin());
vector<RealType>::iterator it(C.begin()),last(C.end());
ValueType psiOld=0.0,psiNew=0.0;
for(int i=0; i<C.size(); i++){
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType tmp2 = C[i]*detValues_up[upC]*detValues_dn[dnC];
psiOld += tmp2;
psiNew += tmp2*detsRatios[dnC];
while(it != last) {
ValueType tmp = (*it)*detValues_up[*upC]*detValues_dn[*dnC];
psiNew += tmp*detsRatios[*dnC];
psiOld += tmp;
it++;upC++;dnC++;
}
curRatio = psiNew/psiOld;
RatioTimer.stop();
RatioTimer.stop();
return curRatio;
}
}
@ -487,6 +499,7 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
{
// this should depend on the type of update, ratio / ratioGrad
// for now is incorrect fot ratio(P,iat,dG,dL) updates
AccRejTimer.start();
if(DetID[iat] == 0) {
for(int i=0; i<dets_up.size(); i++)
dets_up[i]->acceptMove(P,iat);
@ -572,10 +585,12 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
break;
}
}
AccRejTimer.stop();
}
void MultiSlaterDeterminant::restore(int iat)
{
AccRejTimer.start();
if(DetID[iat] == 0) {
for(int i=0; i<dets_up.size(); i++)
dets_up[i]->restore(iat);
@ -584,6 +599,7 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
dets_dn[i]->restore(iat);
}
curRatio=1.0;
AccRejTimer.stop();
}
void MultiSlaterDeterminant::update(ParticleSet& P
@ -787,10 +803,24 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
{
if(Optimizable)
{
for(int i=0; i<C.size(); i++)
{
int loc=myVars.where(i);
if(loc>=0) C[i]=myVars[i]=active[loc];
if(usingCSF) {
for(int i=0; i<CSFcoeff.size()-1; i++) {
int loc=myVars.where(i);
if(loc>=0) CSFcoeff[i+1]=myVars[i]=active[loc];
}
int cnt=0;
for(int i=0; i<DetsPerCSF.size(); i++) {
for(int k=0; k<DetsPerCSF[i]; k++) {
C[cnt] = CSFcoeff[i]*CSFexpansion[cnt];
cnt++;
}
}
} else {
for(int i=0; i<C.size()-1; i++)
{
int loc=myVars.where(i);
if(loc>=0) C[i+1]=myVars[i]=active[loc];
}
}
//for(int i=0; i<SDets.size(); i++) SDets[i]->resetParameters(active);
}
@ -822,43 +852,92 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
// need to modify for CSF later on, right now assume Slater Det basis
if (recalculate)
{
int n = P.getTotalNum();
ValueType psi = std::cos(PhaseValue)*std::exp(LogValue);
ValueType psiinv = 1.0/psi;;
ValueType lapl_sum=0.0;
ParticleSet::ParticleGradient_t g(n);
ValueType gg=0.0, ggP=0.0;
g=0.0;
for(int i=0; i<C.size(); i++){
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType tmp = C[i]*detValues_up[upC]*detValues_dn[dnC]*psiinv;
lapl_sum += tmp*(Sum(lapls_up[upC])+Sum(lapls_dn[dnC]));
g += tmp*grads_up[upC];
g += tmp*grads_dn[dnC];
//cout<<"upC: " <<upC <<" " <<"sum: " <<Sum(lapls_up[upC]) <<" tmp: " <<tmp <<" psiC: " <<psi <<endl;
//cout<<"dnC: " <<upC <<" " <<"sum: " <<Sum(lapls_dn[dnC]) <<endl;
}
gg=Dot(g,g);
ggP=Dot(P.G,g);
if(usingCSF) {
int n = P.getTotalNum();
ValueType psi = std::cos(PhaseValue)*std::exp(LogValue);
ValueType psiinv = 1.0/psi;;
ValueType lapl_sum=0.0;
ParticleSet::ParticleGradient_t g(n);
ValueType gg=0.0, ggP=0.0;
g=0.0;
for(int i=0; i<C.size(); i++){
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType tmp = C[i]*detValues_up[upC]*detValues_dn[dnC]*psiinv;
lapl_sum += tmp*(Sum(lapls_up[upC])+Sum(lapls_dn[dnC]));
g += tmp*grads_up[upC];
g += tmp*grads_dn[dnC];
}
gg=Dot(g,g)-Dot(P.G,g);
// ggP=Dot(P.G,g);
//cout<<"gg: " <<gg <<" ggP:" <<ggP <<endl;
for(int i=0; i<C.size(); i++){
int kk=myVars.where(i);
if (kk<0) continue;
//dlogpsi[kk] = cdet;
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType cdet=detValues_up[upC]*detValues_dn[dnC]*psiinv;
convert(cdet,dlogpsi[kk]);
ValueType dhpsi = (-0.5*cdet)*
( Sum(lapls_up[upC])+Sum(lapls_dn[dnC])-lapl_sum
+2.0*(gg-Dot(g,grads_up[upC])-Dot(g,grads_dn[dnC])
+Dot(P.G,grads_up[upC])+Dot(P.G,grads_dn[dnC])-ggP));
convert(dhpsi,dhpsioverpsi[kk]);
int num=CSFcoeff.size()-1;
int cnt=0;
// this one is not optable
cnt+=DetsPerCSF[0];
int ip(1);
for(int i=0; i<num; i++,ip++) {
int kk=myVars.where(i);
if (kk<0) {
cnt+=DetsPerCSF[ip];
continue;
}
ValueType cdet=0.0,q0=0.0,v1=0.0,v2=0.0;
for(int k=0; k<DetsPerCSF[ip]; k++) {
int upC = C2node_up[cnt];
int dnC = C2node_dn[cnt];
ValueType tmp=CSFexpansion[cnt]*detValues_up[upC]*detValues_dn[dnC]*psiinv;
cdet+=tmp;
q0 += tmp*(Sum(lapls_up[upC])+Sum(lapls_dn[dnC]));
v1 += tmp*(Dot(P.G,grads_up[upC])-Dot(g,grads_up[upC]));
v2 += tmp*(Dot(P.G,grads_dn[dnC])-Dot(g,grads_dn[dnC]));
cnt++;
}
convert(cdet,dlogpsi[kk]);
ValueType dhpsi = -0.5*(q0-cdet*lapl_sum)
-cdet*gg-v1-v2;
//ValueType dhpsi = -0.5*(tmp1*laplSum_up[upC]+tmp2*laplSum_dn[dnC]
// -cdet*lapl_sum)
// -cdet*gg-(tmp1*v1+tmp2*v2);
convert(dhpsi,dhpsioverpsi[kk]);
}
} else {
int n = P.getTotalNum();
ValueType psi = std::cos(PhaseValue)*std::exp(LogValue);
ValueType psiinv = 1.0/psi;;
ValueType lapl_sum=0.0;
ParticleSet::ParticleGradient_t g(n);
ValueType gg=0.0, ggP=0.0;
g=0.0;
for(int i=0; i<C.size(); i++){
int upC = C2node_up[i];
int dnC = C2node_dn[i];
ValueType tmp = C[i]*detValues_up[upC]*detValues_dn[dnC]*psiinv;
lapl_sum += tmp*(Sum(lapls_up[upC])+Sum(lapls_dn[dnC]));
g += tmp*grads_up[upC];
g += tmp*grads_dn[dnC];
}
gg=Dot(g,g);
ggP=Dot(P.G,g);
int ip(1);
for(int i=0; i<C.size()-1; i++,ip++){
int kk=myVars.where(i);
if (kk<0) continue;
//dlogpsi[kk] = cdet;
int upC = C2node_up[ip];
int dnC = C2node_dn[ip];
ValueType cdet=detValues_up[upC]*detValues_dn[dnC]*psiinv;
convert(cdet,dlogpsi[kk]);
ValueType dhpsi = (-0.5*cdet)*
( Sum(lapls_up[upC])+Sum(lapls_dn[dnC])-lapl_sum
+2.0*(gg-Dot(g,grads_up[upC])-Dot(g,grads_dn[dnC])
+Dot(P.G,grads_up[upC])+Dot(P.G,grads_dn[dnC])-ggP));
convert(dhpsi,dhpsioverpsi[kk]);
}
}
}
}
@ -871,7 +950,9 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
RatioAllTimer.reset();
Ratio1AllTimer.reset();
UpdateTimer.reset();
AccRejTimer.reset();
EvaluateTimer.reset();
evalOrbTimer.reset();
TimerManager.addTimer (&RatioTimer);
TimerManager.addTimer (&RatioGradTimer);
TimerManager.addTimer (&RatioAllTimer);
@ -879,7 +960,9 @@ DiracDeterminantBase* adet = new DiracDeterminantBase((SPOSetBasePtr) clone->spo
TimerManager.addTimer (&Ratio1GradTimer);
TimerManager.addTimer (&Ratio1AllTimer);
TimerManager.addTimer (&UpdateTimer);
TimerManager.addTimer (&AccRejTimer);
TimerManager.addTimer (&EvaluateTimer);
TimerManager.addTimer (&evalOrbTimer);
}
}

13
src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h
View File

@ -54,7 +54,7 @@ namespace qmcplusplus
public:
void registerTimers();
NewTimer RatioTimer,RatioGradTimer,RatioAllTimer,UpdateTimer,EvaluateTimer;
NewTimer Ratio1Timer,Ratio1GradTimer,Ratio1AllTimer;
NewTimer Ratio1Timer,Ratio1GradTimer,Ratio1AllTimer,AccRejTimer,evalOrbTimer;
typedef DiracDeterminantBase* DiracDeterminantPtr;
typedef SPOSetBase* SPOSetBasePtr;
@ -166,6 +166,7 @@ namespace qmcplusplus
Vector<ParticleSet::ParticleLaplacian_t> templapl;
ValueType curRatio;
ValueType psiCurrent;
ValueVector_t detsRatios;
ValueVector_t lapl_temp;
GradVector_t grad_temp;
@ -174,6 +175,16 @@ namespace qmcplusplus
ParticleSet::ParticleLaplacian_t myL;
opt_variables_type myVars;
// CSFs
bool usingCSF;
// coefficients of csfs, these are only used during optm
vector<RealType> CSFcoeff;
// number of dets per csf
vector<int> DetsPerCSF;
// coefficient of csf expansion (smaller dimension)
vector<RealType> CSFexpansion;
};
}

26
src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp
View File

@ -27,7 +27,8 @@ namespace qmcplusplus {
Ratio1GradTimer("MultiSlaterDeterminantFast::detEval_ratioGrad"),
Ratio1AllTimer("MultiSlaterDeterminantFast::detEval_ratio(all)"),
UpdateTimer("MultiSlaterDeterminantFast::updateBuffer"),
EvaluateTimer("MultiSlaterDeterminantFast::evaluate")
EvaluateTimer("MultiSlaterDeterminantFast::evaluate"),
AccRejTimer("MultiSlaterDeterminantFast::Accept_Reject")
{
registerTimers();
//Optimizable=true;
@ -310,11 +311,13 @@ namespace qmcplusplus {
vector<RealType>::iterator it(C.begin()),last(C.end());
ValueType psiNew=0.0;
GradType dummy=0.0;
it=C.begin();last=C.end();
while(it != last) {
psiNew += (*it)*detValues_up[*upC]*detValues_dn[*dnC];
dummy += (*it)*grads_up(*upC,iat-N1)*detValues_dn[*dnC];
it++;upC++;dnC++;
}
grad_iat+=dummy/psiNew;
curRatio = psiNew/psiCurrent;
RatioGradTimer.stop();
@ -511,7 +514,7 @@ namespace qmcplusplus {
// for now is incorrect fot ratio(P,iat,dG,dL) updates
// update psiCurrent,myG_temp,myL_temp
AccRejTimer.start();
psiCurrent *= curRatio;
curRatio=1.0;
Dets[DetID[iat]]->acceptMove(P,iat);
@ -525,6 +528,7 @@ namespace qmcplusplus {
default:
break;
}
AccRejTimer.stop();
// Dets[0]->evaluateForWalkerMove(P);
// Dets[1]->evaluateForWalkerMove(P);
@ -566,8 +570,10 @@ namespace qmcplusplus {
void MultiSlaterDeterminantFast::restore(int iat)
{
AccRejTimer.start();
Dets[DetID[iat]]->restore(iat);
curRatio=1.0;
AccRejTimer.stop();
}
void MultiSlaterDeterminantFast::update(ParticleSet& P
@ -620,16 +626,10 @@ namespace qmcplusplus {
OrbitalBase::RealType MultiSlaterDeterminantFast::updateBuffer(ParticleSet& P, BufferType& buf, bool fromscratch)
{
UpdateTimer.start();
if(fromscratch) {
Dets[0]->updateBuffer(P,buf,true);
Dets[1]->updateBuffer(P,buf,true);
} else {
// FIX FIX FIX: right now, I need to allow to recalculate
// dets, grads and lapls without recomputing orbitals,
// for now i always recalculate
Dets[0]->updateBuffer(P,buf,true);
Dets[1]->updateBuffer(P,buf,true);
}
Dets[0]->updateBuffer(P,buf,fromscratch);
Dets[1]->updateBuffer(P,buf,fromscratch);
//Dets[0]->updateBuffer(P,buf,true);
//Dets[1]->updateBuffer(P,buf,true);
// can this change over time??? I don't know yet
ValueVector_t& detValues_up = Dets[0]->detValues;
@ -953,6 +953,7 @@ namespace qmcplusplus {
Ratio1AllTimer.reset();
UpdateTimer.reset();
EvaluateTimer.reset();
AccRejTimer.reset();
TimerManager.addTimer (&RatioTimer);
TimerManager.addTimer (&RatioGradTimer);
TimerManager.addTimer (&RatioAllTimer);
@ -961,6 +962,7 @@ namespace qmcplusplus {
TimerManager.addTimer (&Ratio1AllTimer);
TimerManager.addTimer (&UpdateTimer);
TimerManager.addTimer (&EvaluateTimer);
TimerManager.addTimer (&AccRejTimer);
}

2
src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.h
View File

@ -56,7 +56,7 @@ namespace qmcplusplus
void registerTimers();
NewTimer RatioTimer,RatioGradTimer,RatioAllTimer,UpdateTimer,EvaluateTimer;
NewTimer Ratio1Timer,Ratio1GradTimer,Ratio1AllTimer;
NewTimer Ratio1Timer,Ratio1GradTimer,Ratio1AllTimer, AccRejTimer;
typedef MultiDiracDeterminantBase* DiracDeterminantPtr;
typedef SPOSetBase* SPOSetBasePtr;

257
src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp
View File

@ -595,211 +595,14 @@ namespace qmcplusplus
{
bool success=true;
/*********************************
1. read ci_configurations and coefficients from xml
2. get unique set of determinants
3. create excitation tree for both spin channels
4. build mapping from original expansion to location in the tree
*********************************/
vector<ci_configuration> confgList_up, uniqueConfg_up;
vector<ci_configuration> confgList_dn, uniqueConfg_dn;
ci_configuration baseC_up;
ci_configuration baseC_dn;
vector<RealType>& coeff = multiSD->C;
vector<ci_configuration> uniqueConfg_up, uniqueConfg_dn;
vector<std::string> CItags;
bool optimizeCI;
string optCI="no";
OhmmsAttributeSet ciAttrib;
ciAttrib.add (optCI,"optimize");
ciAttrib.add (optCI,"Optimize");
ciAttrib.put(cur);
bool optimizeCI = optCI=="yes";
xmlNodePtr curRoot=cur,DetListNode;
string cname;
cur = curRoot->children;
while (cur != NULL)//check the basis set
{
getNodeName(cname,cur);
if(cname == "detlist")
{
DetListNode=cur;
app_log() <<"Found determinant list. \n";
}
cur = cur->next;
}
int NCA,NCB,NEA,NEB,nstates,ndets=0,count=0;
string Dettype="DETS";
OhmmsAttributeSet spoAttrib;
spoAttrib.add (NCA, "nca");
spoAttrib.add (NCB, "ncb");
spoAttrib.add (NEA, "nea");
spoAttrib.add (NEB, "neb");
spoAttrib.add (ndets, "size");
spoAttrib.add (Dettype, "type");
spoAttrib.add (nstates, "nstates");
spoAttrib.put(DetListNode);
if(ndets==0) {
APP_ABORT("size==0 in detlist is not allowed. Use slaterdeterminant in this case.\n");
}
if(Dettype != "DETS" && Dettype != "Determinants") {
APP_ABORT("Only allowed type in detlist is DETS. CSF not implemented yet.\n");
}
// cheating until I fix the converter
NCA = multiSD->nels_up-NEA;
NCB = multiSD->nels_dn-NEB;
if(multiSD->nels_up != (NCA+NEA)) {
APP_ABORT("Number of up electrons in ParticleSet doesn't agree with NCA+NEA in detlist.");
}
if(multiSD->nels_dn != (NCB+NEB)) {
APP_ABORT("Number of down electrons in ParticleSet doesn't agree with NCB+NEB in detlist.");
}
if(multiSD->spo_up->refPhi->getOrbitalSetSize() < NCA+nstates) {
APP_ABORT("Number of states in SPOSet is smaller than NCA+nstates in detlist.");
}
if(multiSD->spo_dn->refPhi->getOrbitalSetSize() < NCB+nstates) {
APP_ABORT("Number of states in SPOSet is smaller than NCB+nstates in detlist.");
}
cur = DetListNode->children;
ci_configuration dummyC_alpha;
ci_configuration dummyC_beta;
dummyC_alpha.occup.resize(NCA+nstates,false);
for(int i=0; i<NCA+NEA; i++) dummyC_alpha.occup[i]=true;
dummyC_beta.occup.resize(NCB+nstates,false);
for(int i=0; i<NCB+NEB; i++) dummyC_beta.occup[i]=true;
app_log() <<"alpha reference: \n" <<dummyC_alpha;
app_log() <<"beta reference: \n" <<dummyC_beta;
while (cur != NULL)//check the basis set
{
getNodeName(cname,cur);
if(cname == "configuration" || cname == "ci")
{
RealType ci=0.0;
string alpha,beta,tag;
OhmmsAttributeSet confAttrib;
confAttrib.add(ci,"coeff");
confAttrib.add(alpha,"alpha");
confAttrib.add(beta,"beta");
confAttrib.add(tag,"id");
confAttrib.put(cur);
int nq=0,na,nr;
if(alpha.size() < nstates)
{
cerr<<"alpha: " <<alpha <<endl;
APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates");
}
for(int i=0; i<nstates; i++)
{
if(alpha[i] != '0' && alpha[i] != '1') {
cerr<<alpha <<endl;
APP_ABORT("Found incorrect determinant label.");
}
if(alpha[i] == '1') nq++;
}
if(nq != NEA) {
cerr<<"alpha: " <<alpha <<endl;
APP_ABORT("Found incorrect alpha determinant label. noccup != nca+nea");
}
nq=0;
if(beta.size() < nstates)
{
cerr<<"beta: " <<beta <<endl;
APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates");
}
for(int i=0; i<nstates; i++)
{
if(beta[i] != '0' && beta[i] != '1') {
cerr<<beta <<endl;
APP_ABORT("Found incorrect determinant label.");
}
if(beta[i] == '1') nq++;
}
if(nq != NEB) {
cerr<<"beta: " <<beta <<endl;
APP_ABORT("Found incorrect beta determinant label. noccup != ncb+neb");
}
count++;
coeff.push_back(ci);
CItags.push_back(tag);
confgList_up.push_back(dummyC_alpha);
for(int i=0; i<NCA; i++) confgList_up.back().occup[i]=true;
for(int i=NCA; i<NCA+nstates; i++)
confgList_up.back().occup[i]= (alpha[i-NCA]=='1');
confgList_dn.push_back(dummyC_beta);
for(int i=0; i<NCB; i++) confgList_dn.back().occup[i]=true;
for(int i=NCB; i<NCB+nstates; i++)
confgList_dn.back().occup[i]=(beta[i-NCB]=='1');
}
cur = cur->next;
}
if(count != ndets) {
cerr<<"count, ndets: " <<count <<" " <<ndets <<endl;
APP_ABORT("Problems reading determinant ci_configurations. Found a number of determinants inconsistent with xml file size parameter.\n");
}
if(confgList_up.size() != ndets || confgList_dn.size() != ndets || coeff.size() != ndets) {
APP_ABORT("Problems reading determinant ci_configurations.");
}
multiSD->C2node_up.resize(coeff.size());
multiSD->C2node_dn.resize(coeff.size());
app_log() <<"Found " <<coeff.size() <<" terms in the MSD expansion.\n";
for(int i=0; i<confgList_up.size(); i++)
{
bool found=false;
int k=-1;
for(int j=0; j<uniqueConfg_up.size(); j++)
{
if(confgList_up[i] == uniqueConfg_up[j]) {
found=true;
k=j;
break;
}
}
if(found) {
multiSD->C2node_up[i]=k;
} else {
uniqueConfg_up.push_back(confgList_up[i]);
multiSD->C2node_up[i]=uniqueConfg_up.size()-1;
}
}
for(int i=0; i<confgList_dn.size(); i++)
{
bool found=false;
int k=-1;
for(int j=0; j<uniqueConfg_dn.size(); j++)
{
if(confgList_dn[i] == uniqueConfg_dn[j]) {
found=true;
k=j;
break;
}
}
if(found) {
multiSD->C2node_dn[i]=k;
} else {
uniqueConfg_dn.push_back(confgList_dn[i]);
multiSD->C2node_dn[i]=uniqueConfg_dn.size()-1;
}
}
app_log() <<"Found " <<uniqueConfg_up.size() <<" unique up determinants.\n";
app_log() <<"Found " <<uniqueConfg_dn.size() <<" unique down determinants.\n";
int nels_up = multiSD->nels_up;
int nels_dn = multiSD->nels_dn;
success = readDetList(cur,uniqueConfg_up,uniqueConfg_dn,multiSD->C2node_up, multiSD->C2node_dn,CItags,multiSD->C,optimizeCI,nels_up,nels_dn,multiSD->CSFcoeff,multiSD->DetsPerCSF,multiSD->CSFexpansion,multiSD->usingCSF);
if(!success) return false;
multiSD->resize(uniqueConfg_up.size(),uniqueConfg_dn.size());
SPOSetProxyForMSD* spo = multiSD->spo_up;
@ -833,21 +636,43 @@ namespace qmcplusplus
multiSD->dets_dn.push_back(adet);
}
if (multiSD->CSFcoeff.size()==1 || multiSD->C.size()==1) optimizeCI=false;
if(optimizeCI) {
app_log() <<"CI coefficients are optimizable. \n";
string resetCI("no");
OhmmsAttributeSet spoAttrib;
spoAttrib.add (resetCI, "reset_coeff");
spoAttrib.put(cur);
if (resetCI=="yes")
{
if(multiSD->usingCSF) for(int i=1; i<multiSD->CSFcoeff.size(); i++) multiSD->CSFcoeff[i]=0;
else for(int i=1; i<multiSD->C.size(); i++)multiSD->C[i]=0;
app_log() <<"CI coefficients are reset. \n";
}
multiSD->Optimizable=true;
// multiSD->myVars.insert(CItags[0],coeff[0],false,optimize::LINEAR_P);
for(int i=1; i<coeff.size(); i++) {
//std::stringstream sstr;
//sstr << "CIcoeff" << "_" << i;
multiSD->myVars.insert(CItags[i],coeff[i],true,optimize::LINEAR_P);
if(multiSD->usingCSF) {
// multiSD->myVars.insert(CItags[0],multiSD->CSFcoeff[0],false,optimize::LINEAR_P);
for(int i=1; i<multiSD->CSFcoeff.size(); i++) {
//std::stringstream sstr;
//sstr << "CIcoeff" << "_" << i;
multiSD->myVars.insert(CItags[i],multiSD->CSFcoeff[i],true,optimize::LINEAR_P);
}
} else {
// multiSD->myVars.insert(CItags[0],multiSD->C[0],false,optimize::LINEAR_P);
for(int i=1; i<multiSD->C.size(); i++) {
//std::stringstream sstr;
//sstr << "CIcoeff" << "_" << i;
multiSD->myVars.insert(CItags[i],multiSD->C[i],true,optimize::LINEAR_P);
}
}
}
// //store this for making clones later.
// multiSD->confgList_up=confgList_up;
// multiSD->uniqueConfg_up=uniqueConfg_up;
// multiSD->confgList_dn=confgList_dn;
// multiSD->uniqueConfg_dn=uniqueConfg_dn;
else
{
app_log() <<"CI coefficients are not optimizable. \n";
multiSD->Optimizable=false;
}
return success;
}
@ -930,6 +755,7 @@ namespace qmcplusplus
for(int i=0; i<NCA+NEA; i++) dummyC_alpha.occup[i]=true;
dummyC_beta.occup.resize(NCB+nstates,false);
for(int i=0; i<NCB+NEB; i++) dummyC_beta.occup[i]=true;
RealType sumsq_qc=0.0;
//app_log() <<"alpha reference: \n" <<dummyC_alpha;
//app_log() <<"beta reference: \n" <<dummyC_beta;
@ -950,6 +776,7 @@ namespace qmcplusplus
confAttrib.add(tag,"id");
confAttrib.add(OccString,"occ");
confAttrib.put(cur);
if(qc_ci == 0.0) qc_ci = ci;
if(abs(qc_ci) < cutoff) {
cur = cur->next;
@ -959,6 +786,7 @@ namespace qmcplusplus
cnt0++;
CSFcoeff.push_back(ci);
sumsq_qc += qc_ci*qc_ci;
DetsPerCSF.push_back(0);
CItags.push_back(tag);
count++;
@ -1056,6 +884,7 @@ namespace qmcplusplus
confAttrib.add(beta,"beta");
confAttrib.add(tag,"id");
confAttrib.put(cur);
if(qc_ci == 0.0) qc_ci = ci;
if(abs(qc_ci) < cutoff) {
cur = cur->next;
@ -1105,6 +934,7 @@ namespace qmcplusplus
count++;
coeff.push_back(ci);
sumsq_qc += qc_ci*qc_ci;
CItags.push_back(tag);
confgList_up.push_back(dummyC_alpha);
for(int i=0; i<NCA; i++) confgList_up.back().occup[i]=true;
@ -1135,6 +965,7 @@ namespace qmcplusplus
RealType sumsq=0.0;
for(int i=0; i<coeff.size(); i++) sumsq += coeff[i]*coeff[i];
app_log() <<"Norm of ci vector (sum of ci^2): " <<sumsq <<endl;
app_log() <<"Norm of qchem ci vector (sum of qchem_ci^2): " <<sumsq_qc <<endl;
for(int i=0; i<confgList_up.size(); i++)
{