Commit 251668a2 authored by Chad Hoyer's avatar Chad Hoyer Committed by Hang Hu
Browse files

Added scrFile functionality with -s flag. The scrFile is meant to have data...

Added scrFile functionality with -s flag. The scrFile is meant to have data extracted from it without any connection to rstFile. Save fchk as scrFile with guess=fchk.
parent 1b4e97b3
......@@ -43,9 +43,14 @@ namespace ChronusQ {
// See src/cxxapi/input/*opts.cxx for documentation
// Parse the options relating to the Molecule object
Molecule CQMoleculeOptions(std::ostream &, CQInputFile &);
Molecule CQMoleculeOptions(std::ostream &, CQInputFile &, std::string &);
void CQMOLECULE_VALID(std::ostream&, CQInputFile &);
void parseGeomInp(Molecule &, std::string &, std::ostream &);
void parseGeomFchk(Molecule &, std::string &, std::ostream &);
// Parse the options relating to the BasisSet
std::shared_ptr<BasisSet> CQBasisSetOptions(std::ostream &, CQInputFile &,
Molecule &, std::string);
......
......@@ -228,7 +228,12 @@ namespace ChronusQ {
void RandomGuess();
void ReadGuessMO();
void ReadGuess1PDM();
void FchkGuessMO();
// Fchk-related functions
std::vector<int> fchkToCQMO();
void reorderAngMO(std::vector<int> sl, MatsT* tmo, int sp);
void reorderSpinMO();
// SCF procedural functions (see include/singleslater/scf.hpp for docs)
......
......@@ -46,7 +46,8 @@ namespace ChronusQ {
SAD,
RANDOM,
READMO,
READDEN
READDEN,
FCHKMO
};
/**
......@@ -154,6 +155,9 @@ namespace ChronusQ {
// Save / Restart File
SafeFile savFile;
// Fchk File
std::string fchkFileName;
// Print Controls
size_t printLevel; ///< Print Level
......
This diff is collapsed.
......@@ -160,6 +160,7 @@ namespace ChronusQ {
else if( scfControls.guess == RANDOM ) RandomGuess();
else if( scfControls.guess == READMO ) ReadGuessMO();
else if( scfControls.guess == READDEN ) ReadGuess1PDM();
else if( scfControls.guess == FCHKMO ) FchkGuessMO();
else CErr("Unknown choice for SCF.GUESS",std::cout);
......@@ -820,5 +821,58 @@ namespace ChronusQ {
} // SingleSlater<T>::ReadGuessMO()
/**
* \brief Reads in MOs from fchk file.
*
**/
template <typename MatsT, typename IntsT>
void SingleSlater<MatsT,IntsT>::FchkGuessMO() {
if( printLevel > 0 ){
std::cout << " * Reading in guess orbitals from file " << fchkFileName << "\n";
std::cout << " See documentation of fchkToCQMO() if any problems" << "\n";
}
std::vector<int> shellList;
// Parsing fchk file and populating mo1(mo2)
// Outputs shell listing needed for angular momentum reorganization
shellList=fchkToCQMO();
// Dimension of mo1 and mo2
auto NB = this->nC * basisSet().nBasis;
auto NB2 = NB*NB;
// Reorders shells of mo1 to Chronus ordering
MatsT* mo1tmp = memManager.malloc<MatsT>(NB2);
SetMat('N',NB,NB,MatsT(1.),this->mo[0].pointer(),NB,mo1tmp,NB);
reorderAngMO(shellList,mo1tmp,0);
memManager.free(mo1tmp);
// Reorders shells of mo2 to Chronus ordering
if( this->nC == 1 and not this->iCS){
MatsT* mo2tmp = memManager.malloc<MatsT>(NB2);
SetMat('N',NB,NB,MatsT(1.),this->mo[1].pointer(),NB,mo2tmp,NB);
reorderAngMO(shellList,mo2tmp,1);
memManager.free(mo2tmp);
}
// Reorder spin components
if( this->nC == 2 ) reorderSpinMO();
// Form density from MOs
formDensity();
std::cout << "\n" << std::endl;
if( printLevel > 0 )
std::cout << std::endl
<< " *** Forming Initial Fock Matrix from Guess Density ***\n\n";
std::cout << "\n" << std::endl;
EMPerturbation pert;
formFock(pert,false);
} // SingleSlater<T>::FchkGuessMO()
}; // namespace ChronusQ
......@@ -203,6 +203,7 @@ namespace ChronusQ {
#include <singleslater/extrap.hpp> // Extrapolate header
#include <singleslater/print.hpp> // Print header
#include <singleslater/pop.hpp> // Population analysis
#include <singleslater/fchk.hpp> // Fchk-specific header
#include <singleslater/kohnsham/impl.hpp> // KS headers
#include <singleslater/kohnsham/fxc.hpp> // KS headers
......
......@@ -62,7 +62,7 @@ int main(int argc, char *argv[]) {
} else { // Variable Argc
int c;
while((c = getopt(argc,argv,"i:o:b:z:")) != -1) {
while((c = getopt(argc,argv,"i:o:b:z:s:")) != -1) {
switch(c) {
case('i'):
inFileName = optarg;
......@@ -76,6 +76,9 @@ int main(int argc, char *argv[]) {
case('z'):
oldRstFileName = optarg;
break;
case('s'):
scrFileName = optarg;
break;
default:
abort();
};
......
......@@ -40,6 +40,7 @@ namespace ChronusQ {
std::vector<std::string> allowedKeywords = {
"CHARGE",
"MULT",
"READGEOM",
"GEOM"
};
......@@ -63,7 +64,7 @@ namespace ChronusQ {
*
* \returns Appropriate Molecule object for the input parameters
*/
Molecule CQMoleculeOptions(std::ostream &out, CQInputFile &input) {
Molecule CQMoleculeOptions(std::ostream &out, CQInputFile &input, std::string &scrName) {
Molecule mol;
......@@ -79,14 +80,32 @@ namespace ChronusQ {
//CErr("Unable to set Molecular Spin Multiplicity!", out);
}
// Parse Geometry Read Options
std::string geomReadStr;
try { geomReadStr = input.getData<std::string>("MOLECULE.READGEOM"); }
catch (...) {
geomReadStr="INPUTFILE";
}
// Parse Geometry
std::string geomStr;
try { geomStr = input.getData<std::string>("MOLECULE.GEOM"); }
catch (...) {
CErr("Unable to find Molecular Geometry!", out);
if( geomReadStr == "INPUTFILE" ){
try { geomStr = input.getData<std::string>("MOLECULE.GEOM"); }
catch (...) {
CErr("Unable to find Molecular Geometry in Input File!", out);
}
}
// Parse different files depending on READGEOM
if( geomReadStr == "INPUTFILE" ) parseGeomInp(mol,geomStr,out);
else if( geomReadStr == "FCHK" ) parseGeomFchk(mol,scrName,out);
else CErr("INVALID OPTION FOR READGEOM KEYWORD!");
return mol; // Return Molecule object (no intermediates)
}; // CQMoleculeOptions
void parseGeomInp( Molecule &mol, std::string &geomStr, std::ostream &out ) {
std::istringstream geomStream; geomStream.str(geomStr);
std::vector<std::string> tokens;
......@@ -139,9 +158,125 @@ namespace ChronusQ {
// Output Molecule data
out << mol << std::endl;
return mol; // Return Molecule object (no intermediates)
} // parseGeomInp
}; // CQMoleculeOptions
void parseGeomFchk( Molecule &mol, std::string &fchkName, std::ostream &out ) {
std::ifstream fchkFile;
std::vector<Atom> atoms;
fchkFile.open(fchkName);
if ( fchkFile.good() ){
std::cout << "fchkFile found in parseGeomFchk" << "\n";
}else{
CErr("Could not find fchkFile. Use -s flag.");
}
// Boolean used for fchk reading
bool readAtNum=false, readCoord=false;
// Counter for atomic coordinates
int coordCount=0, atomCount=0;
// TODO: Add non-default isotopes
// Parse fchk file
while( not fchkFile.eof() ) {
std::string line;
std::getline(fchkFile,line);
// Determine position of first and last non-space character
size_t firstNonSpace = line.find_first_not_of(" ");
// Skip blank lines
if( firstNonSpace == std::string::npos ) continue;
// Strip trailing spaces
trim_right(line);
line =
line.substr(firstNonSpace,line.length()-firstNonSpace);
// Split the line into tokens, trim spaces
std::vector<std::string> tokens;
split(tokens,line," ");
for(auto &X : tokens) { trim(X); }
// Finding atomic number in fchk file
if ( tokens.size() > 4 ){
if ( tokens[0] == "Atomic" and tokens[1] == "numbers" ){
readAtNum = true;
continue;
}
}
// Ending atomic number in fchk file
if ( tokens.size() > 4 ){
if ( tokens[0] == "Nuclear" and tokens[1] == "charges" ){
readAtNum = false;
continue;
}
}
// Finding coordinates in fchk file
if ( tokens.size() > 4 ){
if ( tokens[0] == "Current" and tokens[1] == "cartesian" ){
readCoord = true;
continue;
}
}
// Ending atomic number in fchk file
if ( tokens.size() > 5 ){
if ( tokens[2] == "symbols" and tokens[3] == "in" ){
readCoord = false;
continue;
}
}
// Read in atomic numbers
if ( readAtNum ){
for(int i=0; i<tokens.size(); i++){
auto it =
std::find_if(atomicReference.begin(),atomicReference.end(),
[&](std::pair<std::string,Atom> st){
return st.second.atomicNumber == std::stoi(tokens[i]);}
);
// Atomic symbols are currently the first isotope listed in atomicReference
// Trimming everything after the hyphen and passing into defaultIsotope
std::string atmSymb = it->first.substr(0,it->first.find("-",0));
atoms.emplace_back((it == atomicReference.end() ? "X" : defaultIsotope[atmSymb]));
}
} // Read in atomic numbers
// Read in cartesian coordinates from fchk (in Bohr)
if ( readCoord ){
for(int i=0; i<tokens.size(); i++){
atoms[atomCount].coord[coordCount] = std::stod(tokens[i]);
coordCount = coordCount + 1;
// Reset for x, y, z
if( coordCount == 3 ){
coordCount = 0;
atomCount = atomCount + 1;
}
}
} // Read in cartesian coordinates
} // End of fchk file parsing
// Set the Atoms vector in Molecule (calls Molecule::update())
mol.setAtoms(atoms);
// Output Molecule data
out << mol << std::endl;
} // parseGeomFchk
}; // namespace ChronusQ
......@@ -134,6 +134,8 @@ namespace ChronusQ {
ss.scfControls.guess = READMO;
else if( not guessString.compare("READDEN") )
ss.scfControls.guess = READDEN;
else if( not guessString.compare("FCHKMO") )
ss.scfControls.guess = FCHKMO;
else
CErr("Unrecognized entry for SCF.GUESS");
)
......
......@@ -139,7 +139,7 @@ namespace ChronusQ {
// Create Molecule and BasisSet objects
Molecule mol(std::move(CQMoleculeOptions(output,input)));
Molecule mol(std::move(CQMoleculeOptions(output,input,scrFileName)));
std::shared_ptr<BasisSet> basis = CQBasisSetOptions(output,input,mol,"BASIS");
std::shared_ptr<BasisSet> dfbasis = CQBasisSetOptions(output,input,mol,"DFBASIS");
......@@ -158,6 +158,8 @@ namespace ChronusQ {
if( ss->scfControls.guess == READMO or
ss->scfControls.guess == READDEN )
rstExists = true;
if( ss->scfControls.guess == FCHKMO )
ss->fchkFileName = scrFileName;
// Create the restart and scratch files
SafeFile rstFile(rstFileName, rstExists);
......
......@@ -112,9 +112,11 @@ void CONTRACT_TEST(TWOBODY_CONTRACTION_TYPE type, std::string storage) {
// Memory
auto memManager = CQMiscOptions(std::cout,input);
// Dummy scrName since READGEOM=INPUTFILE
std::string scrName;
// Molecule and BasisSet
Molecule mol(std::move(CQMoleculeOptions(std::cout,input)));
Molecule mol(std::move(CQMoleculeOptions(std::cout,input,scrName)));
std::shared_ptr<BasisSet> basis = CQBasisSetOptions(std::cout,input,mol,"BASIS");
// AOIntegrals object
......
......@@ -28,7 +28,7 @@ set(SCF_TEST_BINARY_ROOT "${TEST_BINARY_ROOT}/scf" )
# Set up compilation of SCF test exe
add_executable(scftest ../ut.cxx
rhf.cxx uhf.cxx rohf.cxx x2chf.cxx rhfnr.cxx
rhf.cxx uhf.cxx ghf.cxx rohf.cxx x2chf.cxx rhfnr.cxx
rhf_giao.cxx uhf_giao.cxx rohf_giao.cxx ghf_giao.cxx
rks.cxx uks.cxx x2cks.cxx ks.cxx
misc.cxx ri.cxx
......@@ -49,6 +49,7 @@ file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/rks)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/uhf)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/rohf)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/uks)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/ghf)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/x2c)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/rhf_giao)
file(MAKE_DIRECTORY ${SCF_TEST_BINARY_ROOT}/serial/uhf_giao)
......@@ -84,6 +85,7 @@ include( CQTestGeneration )
# Add tests
add_cq_test( RHF_SCF scftest "RHF.*" )
add_cq_test( UHF_SCF scftest "UHF.*" )
add_cq_test( GHF_SCF scftest "GHF.*" )
add_cq_test( ROHF_SCF scftest "ROHF.*" )
add_cq_test( X2CHF_SCF scftest "X2CHF.*" )
add_cq_test( RHF_NR_SCF scftest "RHF_NR.*" )
......
/*
* This file is part of the Chronus Quantum (ChronusQ) software package
*
* Copyright (C) 2014-2020 Li Research Group (University of Washington)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Contact the Developers:
* E-Mail: xsli@uw.edu
*
*/
#include "scf.hpp"
// KCaKrK sto-3g test for fchk parsing
TEST( GHF, KCaKr_sto3G ) {
CQSCFFCHKTEST( "scf/serial/ghf/KCaKr_sto-3G", "KCaKr_sto-3G.bin.ref", "KCaKr_sto-3G.fchk" );
};
#ifdef _CQ_DO_PARTESTS
#endif
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
......@@ -161,6 +161,130 @@ inline void CQSCFTEST( std::string in, std::string ref,
}
inline void CQSCFFCHKTEST( std::string in, std::string ref, std::string fchk,
double tol = 1e-8,
bool checkSEXP = true,
bool checkSSq = true,
bool checkOctLen = true,
bool checkQuadLen = true,
bool checkDipLen = true,
bool checkEne = true ) {
#ifdef _CQ_GENERATE_TESTS
MPI_Barrier(MPI_COMM_WORLD);
RunChronusQ(TEST_ROOT + in + ".inp","STDOUT",
SCF_TEST_REF + ref, SCF_TEST_REF + fchk);
MPI_Barrier(MPI_COMM_WORLD);
if(MPIRank(MPI_COMM_WORLD) != 0) return;
#else
MPI_Barrier(MPI_COMM_WORLD);
RunChronusQ(TEST_ROOT + in + ".inp","STDOUT",
TEST_OUT + in + ".bin",
SCF_TEST_REF + fchk);
MPI_Barrier(MPI_COMM_WORLD);
if(MPIRank(MPI_COMM_WORLD) != 0) return;
SafeFile refFile(SCF_TEST_REF + ref,true);
SafeFile resFile(TEST_OUT + in + ".bin",true);
double xDummy, yDummy;
std::array<double,3> xDummy3, yDummy3;
std::array<std::array<double,3>,3> xDummy33, yDummy33;
std::array<std::array<std::array<double,3>,3>,3> xDummy333, yDummy333;
/* Check Energy */
if( checkEne ) {
std::cout << " * PERFORMING SCF ENERGY CHECK " << std::endl;
refFile.readData("SCF/TOTAL_ENERGY",&xDummy);
resFile.readData("SCF/TOTAL_ENERGY",&yDummy);
EXPECT_NEAR( xDummy, yDummy, tol ) << "ENERGY TEST FAILED ";
}
/* Check Multipoles */
if( checkDipLen ) {
std::cout << " * PERFORMING SCF DIPOLE (LEN) CHECK " << std::endl;
refFile.readData("SCF/LEN_ELECTRIC_DIPOLE",&xDummy3[0]);
resFile.readData("SCF/LEN_ELECTRIC_DIPOLE",&yDummy3[0]);
for(auto i = 0; i < 3; i++)
EXPECT_NEAR(yDummy3[i], xDummy3[i], tol ) <<
"DIPOLE TEST FAILED IXYZ = " << i;
}
if( checkQuadLen ) {
std::cout << " * PERFORMING SCF QUADRUPOLE (LEN) CHECK " << std::endl;
refFile.readData("SCF/LEN_ELECTRIC_QUADRUPOLE",&xDummy33[0][0]);
resFile.readData("SCF/LEN_ELECTRIC_QUADRUPOLE",&yDummy33[0][0]);
for(auto i = 0; i < 3; i++)
for(auto j = 0; j < 3; j++)
EXPECT_NEAR(yDummy33[i][j], xDummy33[i][j], tol) <<
"QUADRUPOLE TEST FAILED IXYZ = " << i
<< " JXYZ = " << j;
}
if( checkOctLen ) {
std::cout << " * PERFORMING SCF OCTUPOLE (LEN) CHECK " << std::endl;
refFile.readData("SCF/LEN_ELECTRIC_OCTUPOLE",&xDummy333[0][0][0]);
resFile.readData("SCF/LEN_ELECTRIC_OCTUPOLE",&yDummy333[0][0][0]);
for(auto i = 0; i < 3; i++)
for(auto j = 0; j < 3; j++)
for(auto k = 0; k < 3; k++)
EXPECT_NEAR(yDummy333[i][j][k], xDummy333[i][j][k], tol) <<
"OCTUPOLE TEST FAILED IXYZ = " << i
<< " JXYZ = " << j
<< " KXYZ = " << k;
}
/* Check Spin */
if( checkSEXP ) {
std::cout << " * PERFORMING SCF <S> CHECK " << std::endl;
refFile.readData("SCF/S_EXPECT",&xDummy3[0]);
resFile.readData("SCF/S_EXPECT",&yDummy3[0]);
for(auto i = 0; i < 3; i++)
EXPECT_NEAR(yDummy3[i], xDummy3[i], tol ) <<
"<S> TEST FAILED IXYZ = " << i;
}
if( checkSSq ) {
std::cout << " * PERFORMING SCF <S^2> CHECK " << std::endl;
refFile.readData("SCF/S_SQUARED",&xDummy);
resFile.readData("SCF/S_SQUARED",&yDummy);
EXPECT_NEAR(yDummy, xDummy, tol) << "<S^2> TEST FAILED " ;
}
#endif
}
......
[MOLECULE]
charge = 0
mult = 2
# KCaKr
geom:
Kr 0.000000 0.000000 0.000000
K 0.000000 10.000000 0.000000
Ca 0.000000 0.000000 10.000000
[QM]
reference = GHF
job = SCF
[BASIS]
basis = sto-3g
[SCF]
printmos=TRUE
maxiter=2000
guess=fchkmo
[MOLECULE]
charge = 0
mult = 3
# KCaKrK
readgeom=fchk
[QM]
reference = UHF
job = SCF
[BASIS]
basis = sto-3g
[SCF]