Skip to content

Commit

Permalink
fix for spinors in h5 and xml for MSD
Browse files Browse the repository at this point in the history
  • Loading branch information
@camelto2
camelto2 committed Oct 1, 2021
1 parent 3884a95 commit a4ebddb
Show file tree
Hide file tree
Showing 3 changed files with 95 additions and 46 deletions.
65 changes: 37 additions & 28 deletions src/QMCTools/DiracParser.cpp
View file
Expand Up @@ -690,10 +690,8 @@ void DiracParser::parseCOSCIOrbInfo(std::istream& is, const int irrep_idx, OrbTy
}
}

int DiracParser::sortAndStoreCOSCIOrbs(OrbType type)
int DiracParser::sortAndStoreCOSCIOrbs(OrbType type, const int spinor_component)
{
std::vector<int> spinor_component = {0, 2, 1, 3};

int total = 0;
for (int ir = 0; ir < irreps.size(); ir++)
{
Expand All @@ -708,33 +706,30 @@ int DiracParser::sortAndStoreCOSCIOrbs(OrbType type)
}
}

for (int d = 0; d < 4; d++)
for (int i = 0; i < idx.size(); i++)
{
for (int i = 0; i < idx.size(); i++)
if (spinor_component == 0)
{
if (d == 0)
{
EigVal_alpha.push_back(idx[i].first);
total++;
}
int ir = idx[i].second.first;
int mo = idx[i].second.second;
for (int ao = 0; ao < irreps[ir].get_num_ao(); ao++)
EigVec.push_back(irreps[ir].spinor_mo_coeffs[mo][ao][spinor_component[d]]);
EigVal_alpha.push_back(idx[i].first);
total++;
}
//now store KP for this irrep
for (int i = 0; i < idx.size(); i++)
int ir = idx[i].second.first;
int mo = idx[i].second.second;
for (int ao = 0; ao < irreps[ir].get_num_ao(); ao++)
EigVec.push_back(irreps[ir].spinor_mo_coeffs[mo][ao][spinor_component]);
}
//now store KP for this irrep
for (int i = 0; i < idx.size(); i++)
{
if (spinor_component == 0)
{
if (d == 0)
{
EigVal_alpha.push_back(idx[i].first);
total++;
}
int ir = idx[i].second.first;
int mo = idx[i].second.second;
for (int ao = 0; ao < irreps[ir].get_num_ao(); ao++)
EigVec.push_back(kp_irreps[ir].spinor_mo_coeffs[mo][ao][spinor_component[d]]);
EigVal_alpha.push_back(idx[i].first);
total++;
}
int ir = idx[i].second.first;
int mo = idx[i].second.second;
for (int ao = 0; ao < irreps[ir].get_num_ao(); ao++)
EigVec.push_back(kp_irreps[ir].spinor_mo_coeffs[mo][ao][spinor_component]);
}
}

Expand Down Expand Up @@ -786,9 +781,23 @@ void DiracParser::getCOSCI(std::istream& is)
EigVec.clear();
EigVal_alpha.clear();
EigVal_beta.clear();
int total_core = sortAndStoreCOSCIOrbs(OrbType::CORE);
int total_active = sortAndStoreCOSCIOrbs(OrbType::ACTIVE);
int total_virtual = sortAndStoreCOSCIOrbs(OrbType::VIRTUAL);
//save real part of up component of spinor in EigVec first
int total_core = sortAndStoreCOSCIOrbs(OrbType::CORE, 0);
int total_active = sortAndStoreCOSCIOrbs(OrbType::ACTIVE, 0);
int total_virtual = sortAndStoreCOSCIOrbs(OrbType::VIRTUAL, 0);
//save real part of dn compeonent of spinor in Eigvec
sortAndStoreCOSCIOrbs(OrbType::CORE, 2);
sortAndStoreCOSCIOrbs(OrbType::ACTIVE, 2);
sortAndStoreCOSCIOrbs(OrbType::VIRTUAL, 2);
//save imag part of up compeonent of spinor in Eigvec
sortAndStoreCOSCIOrbs(OrbType::CORE, 1);
sortAndStoreCOSCIOrbs(OrbType::ACTIVE, 1);
sortAndStoreCOSCIOrbs(OrbType::VIRTUAL, 1);
//save imag part of dn compeonent of spinor in Eigvec
sortAndStoreCOSCIOrbs(OrbType::CORE, 3);
sortAndStoreCOSCIOrbs(OrbType::ACTIVE, 3);
sortAndStoreCOSCIOrbs(OrbType::VIRTUAL, 3);


//set occstrs for core and virtual
std::string core_occstr;
Expand Down
2 changes: 1 addition & 1 deletion src/QMCTools/DiracParser.h
View file
Expand Up @@ -96,7 +96,7 @@ class DiracParser : public QMCGaussianParserBase, public OhmmsAsciiParser
std::vector<cosciRep> cosciReps;

void parseCOSCIOrbInfo(std::istream& is, const int irrep_idx, OrbType type);
int sortAndStoreCOSCIOrbs(OrbType type);
int sortAndStoreCOSCIOrbs(OrbType type, const int spinor_component);
};


Expand Down
74 changes: 57 additions & 17 deletions src/QMCTools/QMCGaussianParserBase.cpp
View file
Expand Up @@ -931,14 +931,18 @@ void QMCGaussianParserBase::createSPOSetsH5(xmlNodePtr spoUP, xmlNodePtr spoDN)
xmlNewProp(spoUP, (const xmlChar*)"size", (const xmlChar*)nstates_alpha.str().c_str());

//create a spoDN
xmlNewProp(spoDN, (const xmlChar*)"name", (const xmlChar*)"spo-dn");
xmlNewProp(spoDN, (const xmlChar*)"size", (const xmlChar*)nstates_beta.str().c_str());
if (!isSpinor)
{
xmlNewProp(spoDN, (const xmlChar*)"name", (const xmlChar*)"spo-dn");
xmlNewProp(spoDN, (const xmlChar*)"size", (const xmlChar*)nstates_beta.str().c_str());
}


if (DoCusp == true)
{
xmlNewProp(spoUP, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-up.cuspInfo.xml");
xmlNewProp(spoDN, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-dn.cuspInfo.xml");
if (!isSpinor)
xmlNewProp(spoDN, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-dn.cuspInfo.xml");
}


Expand All @@ -962,32 +966,68 @@ void QMCGaussianParserBase::createSPOSetsH5(xmlNodePtr spoUP, xmlNodePtr spoDN)
n++;
}


hout.write(Ctemp, "eigenset_0");


//add occupation DN
occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
xmlAddChild(spoDN, occ_data);

coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
if (SpinRestricted)
xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
else
if (isSpinor)
{
xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"1");
n = numMO * SizeOfBasisSet;
for (int i = 0; i < numMO; i++)
{
for (int j = 0; j < SizeOfBasisSet; j++)
{
Ctemp[i][j] = EigVec[n];
n++;
}
}
hout.write(Ctemp, "eigenset_1");
n = 2 * numMO * SizeOfBasisSet;
for (int i = 0; i < numMO; i++)
{
for (int j = 0; j < SizeOfBasisSet; j++)
{
Ctemp[i][j] = EigVec[n];
n++;
}
}
hout.write(Ctemp, "eigenset_0_imag");
n = 3 * numMO * SizeOfBasisSet;
for (int i = 0; i < numMO; i++)
{
for (int j = 0; j < SizeOfBasisSet; j++)
{
Ctemp[i][j] = EigVec[n];
n++;
}
}
hout.write(Ctemp, "eigenset_1_imag");

hout.write(EigVal_alpha, "eigenval_0");
}
else
{
//add occupation DN
occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
xmlAddChild(spoDN, occ_data);

coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
if (SpinRestricted)
xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
else
{
xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"1");
n = numMO * SizeOfBasisSet;
for (int i = 0; i < numMO; i++)
for (int j = 0; j < SizeOfBasisSet; j++)
{
Ctemp[i][j] = EigVec[n];
n++;
}
hout.write(Ctemp, "eigenset_1");
}
xmlAddChild(spoDN, coeff_data);
}
xmlAddChild(spoDN, coeff_data);

hout.close();
}
Expand Down

0 comments on commit a4ebddb

Please sign in to comment.