QMCPACK
Public Member Functions | List of all members
GaussianFCHKParser Class Reference
+ Inheritance diagram for GaussianFCHKParser:
+ Collaboration diagram for GaussianFCHKParser:

Public Member Functions

 GaussianFCHKParser ()
 
 GaussianFCHKParser (int argc, char **argv)
 
void parse (const std::string &fname) override
 
void getGeometry (std::istream &is)
 
void getGaussianCenters (std::istream &is)
 
- Public Member Functions inherited from QMCGaussianParserBase
 QMCGaussianParserBase ()
 
 QMCGaussianParserBase (int argc, char **argv)
 
virtual ~QMCGaussianParserBase ()=default
 
void setOccupationNumbers ()
 
void createGridNode (int argc, char **argv)
 
void createSPOSets (xmlNodePtr, xmlNodePtr)
 
void createSPOSetsH5 (xmlNodePtr, xmlNodePtr)
 
void PrepareSPOSetsFromH5 (xmlNodePtr, xmlNodePtr)
 
xmlNodePtr createElectronSet (const std::string &ion_tag)
 
xmlNodePtr createIonSet ()
 
xmlNodePtr createCell ()
 
xmlNodePtr createHamiltonian (const std::string &ion_tag, const std::string &psi_tag)
 
xmlNodePtr createBasisSet ()
 
xmlNodePtr createBasisSetWithHDF5 ()
 
xmlNodePtr createCenter (int iat, int _off)
 
void createCenterH5 (int iat, int _off, int numelem)
 
void createShell (int n, int ig, int off_, xmlNodePtr abasis)
 
void createShellH5 (int n, int ig, int off_, int numelem)
 
xmlNodePtr createDeterminantSet ()
 
xmlNodePtr createMultiDeterminantSet ()
 
xmlNodePtr createDeterminantSetWithHDF5 ()
 
xmlNodePtr createMultiDeterminantSetFromH5 ()
 
xmlNodePtr createMultiDeterminantSetCIHDF5 ()
 
xmlNodePtr PrepareDeterminantSetFromHDF5 ()
 
xmlNodePtr createJ3 ()
 
xmlNodePtr createJ2 ()
 
xmlNodePtr createJ1 ()
 
xmlNodePtr parameter (xmlNodePtr Parent, std::string Mypara, std::string a)
 
int numberOfExcitationsCSF (std::string &)
 
virtual void dumpPBC (const std::string &psi_tag, const std::string &ion_tag)
 
virtual void dump (const std::string &psi_tag, const std::string &ion_tag)
 
void dumpStdInput (const std::string &psi_tag, const std::string &ion_tag)
 
void dumpStdInputProd (const std::string &psi_tag, const std::string &ion_tag)
 
- Public Member Functions inherited from OhmmsAsciiParser
void skiplines (std::istream &is, int n)
 
template<class T >
void getValue (std::istream &is, T &aval)
 
template<class T1 , class T2 >
void getValue (std::istream &is, T1 &aval, T2 &bval)
 
template<class IT >
void getValues (std::istream &is, IT first, IT last)
 
int search (std::istream &is, const std::string &keyword)
 
int search (std::istream &is, const std::string &keyword, std::string &the_line)
 
bool lookFor (std::istream &is, const std::string &keyword)
 
bool lookFor (std::istream &is, const std::string &keyword, std::string &the_line)
 

Additional Inherited Members

- Public Types inherited from QMCGaussianParserBase
using value_type = double
 
using SingleParticlePos = ParticleSet::SingleParticlePos
 
- Static Public Member Functions inherited from QMCGaussianParserBase
static void init ()
 
- Public Attributes inherited from QMCGaussianParserBase
bool multideterminant
 
bool multidetH5
 
bool BohrUnit
 
bool SpinRestricted
 
bool Periodicity
 
bool UseHDF5
 
bool PBC
 
bool production
 
bool zeroCI
 
bool orderByExcitation
 
bool addJastrow
 
bool addJastrow3Body
 
bool ECP
 
bool debug
 
bool Structure
 
bool DoCusp
 
bool FixValence
 
bool singledetH5
 
bool optDetCoeffs
 
bool usingCSF
 
bool isSpinor
 
int IonChargeIndex
 
int ValenceChargeIndex
 
int AtomicNumberIndex
 
int NumberOfAtoms
 
int NumberOfEls
 
int target_state
 
int SpinMultiplicity
 
int NumberOfAlpha
 
int NumberOfBeta
 
int SizeOfBasisSet
 
int numMO
 
int readNO
 
int readGuess
 
int numMO2print
 
int ci_size
 
int ci_nca
 
int ci_ncb
 
int ci_nea
 
int ci_neb
 
int ci_nstates
 
int NbKpts
 
int nbexcitedstates
 
double ci_threshold
 
std::vector< double > STwist_Coord
 
std::string Title
 
std::string basisType
 
std::string basisName
 
std::string Normalized
 
std::string CurrentCenter
 
std::string outputFile
 
std::string angular_type
 
std::string expandYlm
 
std::string h5file
 
std::string multih5file
 
std::string WFS_name
 
std::string CodeName
 
const SimulationCell simulation_cell
 
ParticleSet IonSystem
 
std::vector< std::string > GroupName
 
std::vector< int > gShell
 
std::vector< int > gNumber
 
std::vector< int > gBound
 
std::vector< int > Occ_alpha
 
std::vector< int > Occ_beta
 
std::vector< value_type > Qv
 
std::vector< value_type > gExp
 
std::vector< value_type > gC0
 
std::vector< value_type > gC1
 
std::vector< value_type > EigVal_alpha
 
std::vector< value_type > EigVal_beta
 
std::vector< value_type > EigVec
 
std::unique_ptr< xmlNode, void(*)(xmlNodePtr)> gridPtr
 
std::vector< std::string > CIalpha
 
std::vector< std::string > CIbeta
 
std::vector< std::string > CSFocc
 
std::vector< std::vector< std::string > > CSFalpha
 
std::vector< std::vector< std::string > > CSFbeta
 
std::vector< std::vector< double > > CSFexpansion
 
std::vector< double > CIcoeff
 
std::vector< double > X
 
std::vector< double > Y
 
std::vector< double > Z
 
std::vector< int > Image
 
std::vector< int > CIexcitLVL
 
std::vector< std::pair< int, double > > coeff2csf
 
- Public Attributes inherited from OhmmsAsciiParser
char dbuffer [bufferSize]
 
std::vector< std::string > currentWords
 
- Static Public Attributes inherited from QMCGaussianParserBase
static std::map< int, std::string > IonName
 
static std::vector< std::string > gShellType
 
static std::vector< int > gShellID
 
static const std::vector< double > gCoreTable
 
- Static Public Attributes inherited from OhmmsAsciiParser
static const int bufferSize = 200
 

Detailed Description

Definition at line 24 of file GaussianFCHKParser.h.

Constructor & Destructor Documentation

◆ GaussianFCHKParser() [1/2]

GaussianFCHKParser ( )

Definition at line 25 of file GaussianFCHKParser.cpp.

References QMCGaussianParserBase::basisName, QMCGaussianParserBase::FixValence, and QMCGaussianParserBase::Normalized.

26 {
27  basisName = "Gaussian-G2";
28  Normalized = "no";
29  FixValence = true;
30 }

◆ GaussianFCHKParser() [2/2]

GaussianFCHKParser ( int  argc,
char **  argv 
)

Definition at line 32 of file GaussianFCHKParser.cpp.

References QMCGaussianParserBase::basisName, QMCGaussianParserBase::FixValence, and QMCGaussianParserBase::Normalized.

32  : QMCGaussianParserBase(argc, argv)
33 {
34  basisName = "Gaussian-G2";
35  Normalized = "no";
36  FixValence = true;
37 }

Member Function Documentation

◆ getGaussianCenters()

void getGaussianCenters ( std::istream &  is)

Definition at line 608 of file GaussianFCHKParser.cpp.

References QMCGaussianParserBase::gBound, QMCGaussianParserBase::gC0, QMCGaussianParserBase::gC1, OhmmsAsciiParser::getValues(), QMCGaussianParserBase::gExp, QMCGaussianParserBase::gNumber, QMCGaussianParserBase::gShell, qmcplusplus::n, QMCGaussianParserBase::NumberOfAtoms, and OhmmsAsciiParser::search().

Referenced by parse().

609 {
610  //map between Gaussian to Casino Shell notation
611  std::map<int, int> gsMap;
612  gsMap[0] = 1; //s
613  gsMap[-1] = 2; //sp
614  gsMap[1] = 3; //p
615  gsMap[-2] = 4; //d
616  gsMap[-3] = 5; //f
617  gsMap[-4] = 6; //g
618  gsMap[-5] = 7; //l=5 h
619  gsMap[-6] = 8; //l=6 h1??
620  gsMap[-7] = 9; //l=7 h2??
621  gsMap[-8] = 10; //l=8 h3??
622  std::vector<int> n(gShell.size()), dn(NumberOfAtoms, 0);
623  bool SPshell(false);
624  getValues(is, n.begin(), n.end());
625  for (int i = 0; i < n.size(); i++)
626  {
627  gShell[i] = gsMap[n[i]];
628  if (n[i] == -1)
629  SPshell = true;
630  }
631  search(is, "Number");
632  getValues(is, gNumber.begin(), gNumber.end());
633  search(is, "Shell");
634  getValues(is, n.begin(), n.end());
635  for (int i = 0; i < n.size(); i++)
636  dn[n[i] - 1] += 1;
637  gBound[0] = 0;
638  for (int i = 0; i < NumberOfAtoms; i++)
639  {
640  gBound[i + 1] = gBound[i] + dn[i];
641  }
642  search(is, "Primitive");
643  getValues(is, gExp.begin(), gExp.end());
644  search(is, "Contraction");
645  getValues(is, gC0.begin(), gC0.end());
646  if (SPshell)
647  {
648  search(is, "P(S=P)");
649  getValues(is, gC1.begin(), gC1.end());
650  }
651 }
QMCGaussianParserBase::gC0
std::vector< value_type > gC0
Definition: QMCGaussianParserBase.h:103
OhmmsAsciiParser::search
int search(std::istream &is, const std::string &keyword)
Definition: SimpleParser.h:89
QMCGaussianParserBase::gExp
std::vector< value_type > gExp
Definition: QMCGaussianParserBase.h:103
OhmmsAsciiParser::getValues
void getValues(std::istream &is, IT first, IT last)
Definition: SimpleParser.h:76
QMCGaussianParserBase::gNumber
std::vector< int > gNumber
Definition: QMCGaussianParserBase.h:100
QMCGaussianParserBase::gC1
std::vector< value_type > gC1
Definition: QMCGaussianParserBase.h:103

◆ getGeometry()

void getGeometry ( std::istream &  is)

Definition at line 578 of file GaussianFCHKParser.cpp.

References QMCGaussianParserBase::AtomicNumberIndex, ParticleSet::getSpeciesSet(), OhmmsAsciiParser::getValues(), ParticleSet::GroupID, QMCGaussianParserBase::GroupName, QMCGaussianParserBase::IonChargeIndex, QMCGaussianParserBase::IonName, QMCGaussianParserBase::IonSystem, QMCGaussianParserBase::NumberOfAtoms, OHMMS_DIM, ParticleSet::R, and OhmmsAsciiParser::search().

Referenced by parse().

579 {
580  //atomic numbers
581  std::vector<int> atomic_number(NumberOfAtoms);
582  std::vector<double> q(NumberOfAtoms);
583  //read atomic numbers
584  search(is, "Atomic numbers"); //search for Atomic numbers
585  getValues(is, atomic_number.begin(), atomic_number.end());
586  std::streampos pivot = is.tellg();
587  //read effective nuclear charges
588  search(is, "Nuclear"); //search for Nuclear
589  getValues(is, q.begin(), q.end());
590  is.seekg(pivot); //rewind it
591  search(is, "coordinates");
592  std::vector<double> pos(NumberOfAtoms * OHMMS_DIM);
593  getValues(is, pos.begin(), pos.end());
594  SpeciesSet& species(IonSystem.getSpeciesSet());
595  for (int i = 0, ii = 0; i < NumberOfAtoms; i++)
596  {
597  IonSystem.R[i][0] = pos[ii++];
598  IonSystem.R[i][1] = pos[ii++];
599  IonSystem.R[i][2] = pos[ii++];
600  GroupName[i] = IonName[atomic_number[i]];
601  int speciesID = species.addSpecies(GroupName[i]);
602  IonSystem.GroupID[i] = speciesID;
603  species(AtomicNumberIndex, speciesID) = atomic_number[i];
604  species(IonChargeIndex, speciesID) = q[i];
605  }
606 }
QMCGaussianParserBase::GroupName
std::vector< std::string > GroupName
Definition: QMCGaussianParserBase.h:98
OHMMS_DIM
#define OHMMS_DIM
Definition: config.h:64
qmcplusplus::ParticleSet::GroupID
ParticleIndex GroupID
Species ID.
Definition: ParticleSet.h:77
OhmmsAsciiParser::search
int search(std::istream &is, const std::string &keyword)
Definition: SimpleParser.h:89
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
QMCGaussianParserBase::IonName
static std::map< int, std::string > IonName
Definition: QMCGaussianParserBase.h:171
qmcplusplus::ParticleSet::getSpeciesSet
SpeciesSet & getSpeciesSet()
retrun the SpeciesSet of this particle set
Definition: ParticleSet.h:231
OhmmsAsciiParser::getValues
void getValues(std::istream &is, IT first, IT last)
Definition: SimpleParser.h:76
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33

◆ parse()

void parse ( const std::string &  fname)
overridevirtual

Implements QMCGaussianParserBase.

Definition at line 39 of file GaussianFCHKParser.cpp.

References qmcplusplus::abs(), QMCGaussianParserBase::ci_nca, QMCGaussianParserBase::ci_ncb, QMCGaussianParserBase::ci_nea, QMCGaussianParserBase::ci_neb, QMCGaussianParserBase::ci_nstates, QMCGaussianParserBase::ci_threshold, QMCGaussianParserBase::CIalpha, QMCGaussianParserBase::CIbeta, QMCGaussianParserBase::CIcoeff, copy(), ParticleSet::create(), OhmmsAsciiParser::currentWords, BLAS::done, QMCGaussianParserBase::EigVal_alpha, QMCGaussianParserBase::EigVal_beta, QMCGaussianParserBase::EigVec, QMCGaussianParserBase::gBound, QMCGaussianParserBase::gC0, QMCGaussianParserBase::gC1, getGaussianCenters(), getGeometry(), OhmmsAsciiParser::getValues(), getwords(), QMCGaussianParserBase::gExp, QMCGaussianParserBase::gNumber, QMCGaussianParserBase::GroupName, QMCGaussianParserBase::gShell, QMCGaussianParserBase::IonSystem, OhmmsAsciiParser::lookFor(), QMCGaussianParserBase::multideterminant, QMCGaussianParserBase::NumberOfAlpha, QMCGaussianParserBase::NumberOfAtoms, QMCGaussianParserBase::NumberOfBeta, QMCGaussianParserBase::NumberOfEls, QMCGaussianParserBase::numMO, QMCGaussianParserBase::outputFile, parsewords(), OhmmsAsciiParser::search(), QMCGaussianParserBase::SizeOfBasisSet, QMCGaussianParserBase::SpinMultiplicity, QMCGaussianParserBase::SpinRestricted, and QMCGaussianParserBase::Title.

40 {
41  std::ifstream fin(fname.c_str());
42  getwords(currentWords, fin); //1
43  Title = currentWords[0];
44  getwords(currentWords, fin); //2 SP RHF Gen
45  //if(currentWords[1]=="ROHF" || currentWords[1]=="UHF") {
46  if (currentWords[1] == "UHF")
47  {
48  // mmorales: this should be determined by the existence of "Beta MO", since
49  // there are many ways to get unrestricted runs without UHF (e.g. UMP2,UCCD,etc)
50  SpinRestricted = false;
51  //std::cout << " Spin Unrestricted Calculation (UHF). " << std::endl;
52  }
53  else
54  {
55  SpinRestricted = true;
56  //std::cout << " Spin Restricted Calculation (RHF). " << std::endl;
57  }
58  getwords(currentWords, fin); //3 Number of atoms
59  NumberOfAtoms = atoi(currentWords.back().c_str());
60  // TDB: THIS FIX SHOULD BE COMPATIBLE WITH MY OLD FCHK FILES
61  //getwords(currentWords,fin); //4 Charge
62  bool notfound = true;
63  while (notfound)
64  {
65  std::string aline;
66  getwords(currentWords, fin);
67  for (int i = 0; i < currentWords.size(); i++)
68  {
69  if ("Charge" == currentWords[i])
70  {
71  notfound = false;
72  }
73  }
74  }
75  getwords(currentWords, fin); //5 Multiplicity
76  SpinMultiplicity = atoi(currentWords.back().c_str());
77  getwords(currentWords, fin); //6 Number of electrons
78  NumberOfEls = atoi(currentWords.back().c_str());
79  getwords(currentWords, fin); //7 Number of alpha electrons
80  int nup = atoi(currentWords.back().c_str());
81  getwords(currentWords, fin); //8 Number of beta electrons
82  int ndown = atoi(currentWords.back().c_str());
83  //NumberOfAlpha=nup-ndown;
84  NumberOfAlpha = nup;
85  NumberOfBeta = ndown;
86  getwords(currentWords, fin); //9 Number of basis functions
87  SizeOfBasisSet = atoi(currentWords.back().c_str());
88  getwords(currentWords, fin); //10 Number of independent functions
89  int NumOfIndOrb = atoi(currentWords.back().c_str());
90  std::cout << "Number of independent orbitals: " << NumOfIndOrb << std::endl;
91  numMO = NumOfIndOrb;
92  // TDB: THIS ADDITION SHOULD BE COMPATIBLE WITH MY OLD FCHK FILES
93  std::streampos pivottdb = fin.tellg();
94  int ng = 0;
95  notfound = true;
96  while (notfound)
97  {
98  std::string aline;
99  getwords(currentWords, fin);
100  for (int i = 0; i < currentWords.size(); i++)
101  {
102  if ("contracted" == currentWords[i])
103  {
104  ng = atoi(currentWords.back().c_str());
105  notfound = false;
106  }
107  }
108  }
109  // TDB: THIS FIX SHOULD BE COMPATIBLE WITH MY OLD FCHK FILES
110  // getwords(currentWords,fin); //Number of contracted shells
111  // getwords(currentWords,fin); //Number of contracted shells
112  // getwords(currentWords,fin); //Number of contracted shells
113  // int nx=atoi(currentWords.back().c_str()); //Number of exponents
114  int nx = 0;
115  notfound = true;
116  while (notfound)
117  {
118  std::string aline;
119  getwords(currentWords, fin);
120  for (int i = 0; i < currentWords.size(); i++)
121  {
122  if ("primitive" == currentWords[i])
123  {
124  nx = atoi(currentWords.back().c_str()); //Number of exponents
125  notfound = false;
126  }
127  }
128  }
129  //allocate everything here
130  IonSystem.create({NumberOfAtoms});
131  GroupName.resize(NumberOfAtoms);
132  gBound.resize(NumberOfAtoms + 1);
133  gShell.resize(ng);
134  gNumber.resize(ng);
135  gExp.resize(nx);
136  gC0.resize(nx);
137  gC1.resize(nx);
138  // TDB: THIS ADDITION SHOULD BE COMPATIBLE WITH MY OLD FCHK FILES
139  fin.seekg(pivottdb); //rewind it
140  getGeometry(fin);
141  std::cout << "Number of gaussians " << ng << std::endl;
142  std::cout << "Number of primitives " << nx << std::endl;
143  std::cout << "Number of atoms " << NumberOfAtoms << std::endl;
144  search(fin, "Shell types");
145  getGaussianCenters(fin);
146  std::cout << " Shell types reading: OK" << std::endl;
147  // mmorales: SizeOfBasisSet > numMO always, so leave it like this
148  EigVal_alpha.resize(SizeOfBasisSet);
149  EigVal_beta.resize(SizeOfBasisSet);
150  // mmorales HACK HACK HACK, look for a way to rewind w/o closing/opening a file
151  SpinRestricted = !(lookFor(fin, "Beta MO"));
152  fin.close();
153  fin.open(fname.c_str());
154  search(fin, "Alpha Orbital"); //search "Alpha Orbital Energies"
155  // only read NumOfIndOrb
156  getValues(fin, EigVal_alpha.begin(), EigVal_alpha.begin() + numMO);
157  std::cout << " Orbital energies reading: OK" << std::endl;
158  if (SpinRestricted)
159  {
160  EigVec.resize(2 * numMO * SizeOfBasisSet);
161  EigVal_beta = EigVal_alpha;
162  std::vector<value_type> etemp;
163  search(fin, "Alpha MO");
164  getValues(fin, EigVec.begin(), EigVec.begin() + SizeOfBasisSet * numMO);
165  copy(EigVec.begin(), EigVec.begin() + SizeOfBasisSet * numMO, EigVec.begin() + SizeOfBasisSet * numMO);
166  std::cout << " Orbital coefficients reading: OK" << std::endl;
167  }
168  else
169  {
170  EigVec.resize(2 * numMO * SizeOfBasisSet);
171  std::vector<value_type> etemp;
172  search(fin, "Beta Orbital");
173  getValues(fin, EigVal_beta.begin(), EigVal_beta.begin() + numMO);
174  std::cout << " Read Beta Orbital energies: OK" << std::endl;
175  search(fin, "Alpha MO");
176  getValues(fin, EigVec.begin(), EigVec.begin() + SizeOfBasisSet * numMO);
177  search(fin, "Beta MO");
178  getValues(fin, EigVec.begin() + numMO * SizeOfBasisSet, EigVec.begin() + 2 * numMO * SizeOfBasisSet);
179  std::cout << " Alpha and Beta Orbital coefficients reading: OK" << std::endl;
180  }
181  if (multideterminant)
182  {
183  std::ifstream ofile(outputFile.c_str());
184  if (ofile.fail())
185  {
186  std::cerr << "Failed to open output file from gaussian. \n";
187  exit(401);
188  }
189  int statesType = 0;
190  std::streampos beginpos = ofile.tellg();
191  bool found = lookFor(ofile, "SLATER DETERMINANT BASIS");
192  if (!found)
193  {
194  ofile.close();
195  ofile.open(outputFile.c_str());
196  ofile.seekg(beginpos);
197  found = lookFor(ofile, "Slater Determinants");
198  if (found)
199  {
200  statesType = 1;
201  ofile.close();
202  ofile.open(outputFile.c_str());
203  ofile.seekg(beginpos);
204  search(ofile, "Total number of active electrons");
205  getwords(currentWords, ofile);
206  ci_nstates = atoi(currentWords[5].c_str());
207  search(ofile, "Slater Determinants");
208  }
209  else
210  {
211  std::cerr << "Gaussian parser currently works for slater determinant basis only. Use SlaterDet in CAS() or "
212  "improve parser.\n";
213  abort();
214  }
215  }
216  beginpos = ofile.tellg();
217  // mmorales: gaussian by default prints only 50 states
218  // if you want more you need to use, iop(5/33=1), but this prints
219  // a lot of things. So far, I don't know how to change this
220  // without modifying the source code, l510.F or utilam.F
221  found = lookFor(ofile, "Do an extra-iteration for final printing");
222  std::map<int, int> coeff2confg;
223  // 290 FORMAT(1X,7('(',I5,')',F10.7,1X)/(1X,7('(',I5,')',F10.7,1X)))
224  if (found)
225  {
226  // this is tricky, because output depends on diagonalizer used (Davidson or Lanczos) for now hope this works
227  std::streampos tmppos = ofile.tellg();
228  // long output with iop(5/33=1)
229  int coeffType = 0;
230  if (lookFor(ofile, "EIGENVECTOR USED TO COMPUTE DENSITY MATRICES SET"))
231  {
232  coeffType = 1;
233  }
234  else
235  {
236  // short list (50 states) with either Dav or Lanc
237  ofile.close();
238  ofile.open(outputFile.c_str());
239  ofile.seekg(tmppos); //rewind it
240  if (!lookFor(ofile, "EIGENVALUE "))
241  {
242  ofile.close();
243  ofile.open(outputFile.c_str());
244  ofile.seekg(tmppos); //rewind it
245  if (!lookFor(ofile, "Eigenvalue"))
246  {
247  std::cerr << "Failed to find CI voefficients.\n";
248  abort();
249  }
250  }
251  }
252  std::streampos strpos = ofile.tellg();
253  ofile.close();
254  ofile.open(outputFile.c_str());
255  ofile.seekg(beginpos); //rewind it
256  std::string aline;
257  int numCI;
258  if (lookFor(ofile, "NO OF BASIS FUNCTIONS", aline))
259  {
260  parsewords(aline.c_str(), currentWords);
261  numCI = atoi(currentWords[5].c_str());
262  }
263  else
264  {
265  ofile.close();
266  ofile.open(outputFile.c_str());
267  if (lookFor(ofile, "Number of configurations", aline))
268  {
269  parsewords(aline.c_str(), currentWords);
270  numCI = atoi(currentWords[3].c_str());
271  }
272  else
273  {
274  std::cerr << "Problems finding total number of configurations. \n";
275  abort();
276  }
277  }
278  std::cout << "Total number of CI configurations in file (w/o truncation): " << numCI << std::endl;
279  ofile.close();
280  ofile.open(outputFile.c_str());
281  ofile.seekg(strpos); //rewind it
282  int cnt = 0;
283  CIcoeff.clear();
284  CIalpha.clear();
285  CIbeta.clear();
286  if (coeffType == 0)
287  {
288  // short list
289  for (int i = 0; i < 7; i++)
290  {
291  int pos = 2;
292  std::string aline;
293  getline(ofile, aline, '\n');
294  //cout<<"aline: " <<aline << std::endl;
295  for (int i = 0; i < 7; i++)
296  {
297  //int q = atoi( (aline.substr(pos,pos+4)).c_str() );
298  //coeff2confg[q] = cnt++;
299  //CIcoeff.push_back( atof( (aline.substr(pos+6,pos+15)).c_str() ) );
300  //pos+=18;
301  int q = atoi((aline.substr(pos, pos + 7)).c_str());
302  coeff2confg[q] = cnt++;
303  CIcoeff.push_back(atof((aline.substr(pos + 9, pos + 18)).c_str()));
304  pos += 21;
305  //cout<<"confg, coeff: " <<q <<" " <<CIcoeff.back() << std::endl;
306  }
307  }
308  {
309  int pos = 2;
310  std::string aline;
311  getline(ofile, aline, '\n');
312  //int q = atoi( (aline.substr(pos,pos+4)).c_str() );
313  //coeff2confg[q] = cnt++;
314  //CIcoeff.push_back( atof( (aline.substr(pos+6,pos+15)).c_str() ) );
315  int q = atoi((aline.substr(pos, pos + 7)).c_str());
316  coeff2confg[q] = cnt++;
317  CIcoeff.push_back(atof((aline.substr(pos + 9, pos + 18)).c_str()));
318  //cout<<"confg, coeff: " <<q <<" " <<CIcoeff.back() << std::endl;
319  }
320  }
321  else
322  {
323  // long list
324  int nrows = numCI / 5;
325  int nextra = numCI - nrows * 5;
326  int indx[5];
327  ofile.close();
328  ofile.open(outputFile.c_str());
329  ofile.seekg(strpos); //rewind it
330  for (int i = 0; i < nrows; i++)
331  {
332  getwords(currentWords, ofile);
333  if (currentWords.size() != 5)
334  {
335  std::cerr << "Error reading CI configurations-line: " << i << "\n";
336  abort();
337  }
338  for (int k = 0; k < 5; k++)
339  indx[k] = atoi(currentWords[k].c_str());
340  getwords(currentWords, ofile);
341  if (currentWords.size() != 6)
342  {
343  std::cerr << "Error reading CI configurations-line: " << i << "\n";
344  abort();
345  }
346  for (int k = 0; k < 5; k++)
347  {
348  // HACK HACK HACK - how is this done formally???
349  for (int j = 0; j < currentWords[k + 1].size(); j++)
350  if (currentWords[k + 1][j] == 'D')
351  currentWords[k + 1][j] = 'E';
352  double ci = atof(currentWords[k + 1].c_str());
353  if (std::abs(ci) > ci_threshold)
354  {
355  coeff2confg[indx[k]] = cnt++;
356  CIcoeff.push_back(ci);
357  //cout<<"ind,cnt,c: " <<indx[k] <<" " <<cnt-1 <<" " <<CIcoeff.back() << std::endl;
358  }
359  }
360  }
361  if (nextra > 0)
362  {
363  getwords(currentWords, ofile);
364  if (currentWords.size() != nextra)
365  {
366  std::cerr << "Error reading CI configurations last line \n";
367  abort();
368  }
369  for (int k = 0; k < nextra; k++)
370  indx[k] = atoi(currentWords[k].c_str());
371  getwords(currentWords, ofile);
372  if (currentWords.size() != nextra + 1)
373  {
374  std::cerr << "Error reading CI configurations last line \n";
375  abort();
376  }
377  for (int k = 0; k < nextra; k++)
378  {
379  double ci = atof(currentWords[k + 1].c_str());
380  if (std::abs(ci) > ci_threshold)
381  {
382  coeff2confg[indx[k]] = cnt++;
383  CIcoeff.push_back(ci);
384  }
385  }
386  }
387  }
388  std::cout << "Found " << CIcoeff.size() << " coeffficients after truncation. \n";
389  if (statesType == 0)
390  {
391  ofile.close();
392  ofile.open(outputFile.c_str());
393  ofile.seekg(beginpos); //rewind it
394  // this might not work, look for better entry point later
395  search(ofile, "Truncation Level=");
396  //cout<<"found Truncation Level=" << std::endl;
397  getwords(currentWords, ofile);
398  while (!ofile.eof() &&
399  (currentWords[0] != "no." || currentWords[1] != "active" || currentWords[2] != "orbitals"))
400  {
401  // std::cout <<"1. " <<currentWords[0] << std::endl;
402  getwords(currentWords, ofile);
403  }
404  ci_nstates = atoi(currentWords[4].c_str());
405  getwords(currentWords, ofile);
406  // can choose specific irreps if I want...
407  while (currentWords[0] != "Configuration" || currentWords[2] != "Symmetry")
408  {
409  // std::cout <<"2. " <<currentWords[0] << std::endl;
410  getwords(currentWords, ofile);
411  }
412  CIbeta.resize(CIcoeff.size());
413  CIalpha.resize(CIcoeff.size());
414  bool done = false;
415  // can choose specific irreps if I want...
416  while (currentWords[0] == "Configuration" && currentWords[2] == "Symmetry")
417  {
418  int pos = atoi(currentWords[1].c_str());
419  std::map<int, int>::iterator it = coeff2confg.find(pos);
420  //cout<<"3. configuration: " <<currentWords[1].c_str() << std::endl;
421  if (it != coeff2confg.end())
422  {
423  std::string alp(currentWords[4]);
424  std::string beta(currentWords[4]);
425  if (alp.size() != ci_nstates)
426  {
427  std::cerr << "Problem with ci string. \n";
428  abort();
429  }
430  for (int i = 0; i < alp.size(); i++)
431  {
432  if (alp[i] == 'a')
433  alp[i] = '1';
434  if (alp[i] == 'b')
435  alp[i] = '0';
436  if (beta[i] == 'a')
437  beta[i] = '0';
438  if (beta[i] == 'b')
439  beta[i] = '1';
440  }
441  if (done)
442  {
443  // check number of alpha/beta electrons
444  int n1 = 0;
445  for (int i = 0; i < alp.size(); i++)
446  if (alp[i] == '1')
447  n1++;
448  if (n1 != ci_nea)
449  {
450  std::cerr << "Problems with alpha ci string: " << std::endl
451  << alp << std::endl
452  << currentWords[3] << std::endl;
453  abort();
454  }
455  n1 = 0;
456  for (int i = 0; i < beta.size(); i++)
457  if (beta[i] == '1')
458  n1++;
459  if (n1 != ci_neb)
460  {
461  std::cerr << "Problems with beta ci string: " << std::endl
462  << beta << std::endl
463  << currentWords[3] << std::endl;
464  abort();
465  }
466  }
467  else
468  {
469  // count number of alpha/beta electrons
470  ci_nea = 0;
471  for (int i = 0; i < alp.size(); i++)
472  if (alp[i] == '1')
473  ci_nea++;
474  ci_neb = 0;
475  for (int i = 0; i < beta.size(); i++)
476  if (beta[i] == '1')
477  ci_neb++;
478  ci_nca = nup - ci_nea;
479  ci_ncb = ndown - ci_neb;
480  }
481  CIalpha[it->second] = alp;
482  CIbeta[it->second] = beta;
483  //cout<<"alpha: " <<alp <<" - " <<CIalpha[it->second] << std::endl;
484  }
485  getwords(currentWords, ofile);
486  }
487  //cout.flush();
488  }
489  else
490  {
491  // coefficient list obtained with iop(4/21=110)
492  ofile.close();
493  ofile.open(outputFile.c_str());
494  ofile.seekg(beginpos); //rewind it
495  bool first_alpha = true;
496  bool first_beta = true;
497  std::string aline;
498  getline(ofile, aline, '\n');
499  CIbeta.resize(CIcoeff.size());
500  CIalpha.resize(CIcoeff.size());
501  for (int nst = 0; nst < numCI; nst++)
502  {
503  getwords(currentWords, ofile); // state number
504  int pos = atoi(currentWords[0].c_str());
505  std::map<int, int>::iterator it = coeff2confg.find(pos);
506  if (it != coeff2confg.end())
507  {
508  getwords(currentWords, ofile); // state number
509  if (first_alpha)
510  {
511  first_alpha = false;
512  ci_nea = currentWords.size();
513  ci_nca = nup - ci_nea;
514  std::cout << "nca, nea, nstates: " << ci_nca << " " << ci_nea << " " << ci_nstates << std::endl;
515  }
516  else
517  {
518  if (currentWords.size() != ci_nea)
519  {
520  std::cerr << "Problems with alpha string: " << pos << std::endl;
521  abort();
522  }
523  }
524  std::string alp(ci_nstates, '0');
525  for (int i = 0; i < currentWords.size(); i++)
526  {
527  // std::cout <<"i, alpOcc: " <<i <<" " <<atoi(currentWords[i].c_str())-1 << std::endl; std::cout.flush();
528  alp[atoi(currentWords[i].c_str()) - 1] = '1';
529  }
530  getwords(currentWords, ofile); // state number
531  if (first_beta)
532  {
533  first_beta = false;
534  ci_neb = currentWords.size();
535  ci_ncb = ndown - ci_neb;
536  std::cout << "ncb, neb, nstates: " << ci_ncb << " " << ci_neb << " " << ci_nstates << std::endl;
537  }
538  else
539  {
540  if (currentWords.size() != ci_neb)
541  {
542  std::cerr << "Problems with beta string: " << pos << std::endl;
543  abort();
544  }
545  }
546  std::string beta(ci_nstates, '0');
547  for (int i = 0; i < currentWords.size(); i++)
548  {
549  // std::cout <<"i, alpOcc: " <<i <<" " <<atoi(currentWords[i].c_str())-1 << std::endl; std::cout.flush();
550  beta[atoi(currentWords[i].c_str()) - 1] = '1';
551  }
552  CIalpha[it->second] = alp;
553  CIbeta[it->second] = beta;
554  //cout<<"alpha: " <<alp << std::endl <<"beta: " <<beta << std::endl;
555  std::string aline1;
556  getline(ofile, aline1, '\n');
557  }
558  else
559  {
560  getwords(currentWords, ofile);
561  getwords(currentWords, ofile);
562  std::string aline1;
563  getline(ofile, aline1, '\n');
564  }
565  }
566  }
567  ofile.close();
568  }
569  else
570  {
571  std::cerr << "Could not find CI coefficients in gaussian output file. \n";
572  abort();
573  }
574  std::cout << " size of CIalpha,CIbeta: " << CIalpha.size() << " " << CIbeta.size() << std::endl;
575  }
576 }
QMCGaussianParserBase::EigVal_beta
std::vector< value_type > EigVal_beta
Definition: QMCGaussianParserBase.h:104
GaussianFCHKParser::getGaussianCenters
void getGaussianCenters(std::istream &is)
Definition: GaussianFCHKParser.cpp:608
OhmmsAsciiParser::currentWords
std::vector< std::string > currentWords
Definition: SimpleParser.h:49
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
QMCGaussianParserBase::CIcoeff
std::vector< double > CIcoeff
Definition: QMCGaussianParserBase.h:114
QMCGaussianParserBase::GroupName
std::vector< std::string > GroupName
Definition: QMCGaussianParserBase.h:98
copy
void copy(const Array< T1, 3 > &src, Array< T2, 3 > &dest)
Definition: Blitz.h:639
QMCGaussianParserBase::EigVal_alpha
std::vector< value_type > EigVal_alpha
Definition: QMCGaussianParserBase.h:104
QMCGaussianParserBase::gC0
std::vector< value_type > gC0
Definition: QMCGaussianParserBase.h:103
parsewords
unsigned parsewords(const char *inbuf, std::vector< std::string > &slist, const std::string &extra_tokens)
Definition: SimpleParser.cpp:66
QMCGaussianParserBase::EigVec
std::vector< value_type > EigVec
Definition: QMCGaussianParserBase.h:105
OhmmsAsciiParser::search
int search(std::istream &is, const std::string &keyword)
Definition: SimpleParser.h:89
QMCGaussianParserBase::gExp
std::vector< value_type > gExp
Definition: QMCGaussianParserBase.h:103
QMCGaussianParserBase::CIalpha
std::vector< std::string > CIalpha
Definition: QMCGaussianParserBase.h:110
qmcplusplus::ParticleSet::create
void create(const std::vector< int > &agroup)
create grouped particles
Definition: ParticleSet.cpp:125
OhmmsAsciiParser::lookFor
bool lookFor(std::istream &is, const std::string &keyword)
Definition: SimpleParser.h:130
BLAS::done
constexpr double done
Definition: BLAS.hpp:48
OhmmsAsciiParser::getValues
void getValues(std::istream &is, IT first, IT last)
Definition: SimpleParser.h:76
QMCGaussianParserBase::gNumber
std::vector< int > gNumber
Definition: QMCGaussianParserBase.h:100
getwords
int getwords(std::vector< std::string > &slist, std::istream &fp, std::string &aline)
Definition: SimpleParser.cpp:100
QMCGaussianParserBase::CIbeta
std::vector< std::string > CIbeta
Definition: QMCGaussianParserBase.h:110
QMCGaussianParserBase::gC1
std::vector< value_type > gC1
Definition: QMCGaussianParserBase.h:103
GaussianFCHKParser::getGeometry
void getGeometry(std::istream &is)
Definition: GaussianFCHKParser.cpp:578
The documentation for this class was generated from the following files: