Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ROHF guess fix #570

Merged
merged 3 commits into from Jan 15, 2017
Merged

ROHF guess fix #570

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5a069a6
HF: Removes old orbital read/writing.
dgasmith Jan 13, 2017
5b0aba2
HF: Fixes ROHF read intial guess
dgasmith Jan 13, 2017
42f82b6
Merge branch 'master' into rohf_guess
andysim Jan 15, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
276 changes: 8 additions & 268 deletions psi4/src/psi4/libscf_solver/hf.cc
View file
Open in desktop
Expand Up @@ -1388,10 +1388,11 @@ void HF::guess()
nbetapi_ = guess_Cb_->colspi();
nalpha_ = nalphapi_.sum();
nbeta_ = nbetapi_.sum();
doccpi_ = nalphapi_;
soccpi_ = nalphapi_ - nbetapi_;
doccpi_ = nalphapi_ - soccpi_;
}

format_guess();
form_D();

// This is a guess iteration similar to SAD
Expand Down Expand Up @@ -1463,273 +1464,10 @@ void HF::guess()
E_ = 0.0; // don't use this guess in our convergence checks
}

// void HF::save_orbitals()
// {
// psio_->open(PSIF_SCF_MOS,PSIO_OPEN_NEW);

// if (print_)
// outfile->Printf("\n Saving occupied orbitals to File %d.\n", PSIF_SCF_MOS);

// psio_->write_entry(PSIF_SCF_MOS,"SCF ENERGY",(char *) &(E_),sizeof(double));
// psio_->write_entry(PSIF_SCF_MOS,"NIRREP",(char *) &(nirrep_),sizeof(int));
// psio_->write_entry(PSIF_SCF_MOS,"NSOPI",(char *) &(nsopi_[0]),nirrep_*sizeof(int));
// psio_->write_entry(PSIF_SCF_MOS,"NALPHAPI",(char *) &(nalphapi_[0]),nirrep_*sizeof(int));
// psio_->write_entry(PSIF_SCF_MOS,"NBETAPI",(char *) &(nbetapi_[0]),nirrep_*sizeof(int));

// char *basisname = strdup(options_.get_str("BASIS").c_str());
// int basislength = strlen(options_.get_str("BASIS").c_str()) + 1;
// int nbf = basisset_->nbf();

// psio_->write_entry(PSIF_SCF_MOS,"BASIS NAME LENGTH",(char *)(&basislength),sizeof(int));
// psio_->write_entry(PSIF_SCF_MOS,"BASIS NAME",basisname,basislength*sizeof(char));
// psio_->write_entry(PSIF_SCF_MOS,"NUMBER OF BASIS FUNCTIONS",(char *)(&nbf),sizeof(int));

// // upon loading, need to know what value of puream was used
// int old_puream = (basisset_->has_puream() ? 1 : 0);
// psio_->write_entry(PSIF_SCF_MOS,"PUREAM",(char *)(&old_puream),sizeof(int));

// SharedMatrix Ctemp_a(new Matrix("ALPHA MOS", nirrep_, nsopi_, nalphapi_));
// for (int h = 0; h < nirrep_; h++)
// for (int m = 0; m<nsopi_[h]; m++)
// for (int i = 0; i<nalphapi_[h]; i++)
// Ctemp_a->set(h,m,i,Ca_->get(h,m,i));
// Ctemp_a->save(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);

// SharedMatrix Ctemp_b(new Matrix("BETA MOS", nirrep_, nsopi_, nbetapi_));
// for (int h = 0; h < nirrep_; h++)
// for (int m = 0; m<nsopi_[h]; m++)
// for (int i = 0; i<nbetapi_[h]; i++)
// Ctemp_b->set(h,m,i,Cb_->get(h,m,i));
// Ctemp_b->save(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// // Write Fock matrix to file 280 after removing symmetry
// MintsHelper helper(basisset_, options_, 0);
// SharedMatrix sotoao = helper.petite_list()->sotoao();
// SharedMatrix Fa(new Matrix(nbf,nbf));
// SharedMatrix Fb(new Matrix(nbf,nbf));
// Fa->remove_symmetry(Fa_,sotoao);
// Fb->remove_symmetry(Fb_,sotoao);
// Fa->set_name("ALPHA FOCK C1");
// Fb->set_name("BETA FOCK C1");
// Fa->save(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// Fb->save(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);

// psio_->close(PSIF_SCF_MOS,1);
// free(basisname);
// }

// bool HF::do_use_fock_guess() // only use this approach when symmetry changes are causing issues
// {
// psio_->open(PSIF_SCF_MOS,PSIO_OPEN_OLD);
// // Compare current basis with old basis to see if fock guess will work
// bool is_same_basis = false;
// int nbf = basisset_->nbf(), basislength, old_nbf;
// psio_->read_entry(PSIF_SCF_MOS,"BASIS NAME LENGTH",(char *)(&basislength),sizeof(int));
// char *basisnamec = new char[basislength];
// psio_->read_entry(PSIF_SCF_MOS,"BASIS NAME",basisnamec,basislength*sizeof(char));
// psio_->read_entry(PSIF_SCF_MOS,"NUMBER OF BASIS FUNCTIONS",(char *)(&old_nbf),sizeof(int));
// std::string old_basisname(basisnamec); delete[] basisnamec;
// is_same_basis = (options_.get_str("BASIS") == old_basisname && nbf == old_nbf);
// // Compare number of irreps with old number to see if fock guess is necessary
// bool is_different_symmetry = false;
// int old_nirrep;
// psio_->read_entry(PSIF_SCF_MOS,"NIRREP",(char *) &(old_nirrep),sizeof(int));
// is_different_symmetry = (nirrep_ != old_nirrep);
// psio_->close(PSIF_SCF_MOS,1);

// // Final comparison: Use a Fock guess if you are using the same basis as before and
// // the symmetry has changed
// return (is_same_basis && is_different_symmetry);
// }

// void HF::load_fock()
// {
// int nbf = basisset_->nbf();
// psio_->open(PSIF_SCF_MOS,PSIO_OPEN_OLD);
// // Read Fock matrix from file 280, applying current symmetry
// MintsHelper helper(basisset_, options_, 0);
// SharedMatrix aotoso = helper.petite_list()->aotoso();
// SharedMatrix Fa(new Matrix("ALPHA FOCK C1",nbf,nbf));
// SharedMatrix Fb(new Matrix("BETA FOCK C1",nbf,nbf));
// Fa->load(psio_,PSIF_SCF_MOS,Matrix::SubBlocks);
// Fb->load(psio_,PSIF_SCF_MOS,Matrix::SubBlocks);
// Fa_->apply_symmetry(Fa,aotoso);
// Fb_->apply_symmetry(Fb,aotoso);
// psio_->close(PSIF_SCF_MOS,1);
// }

// void HF::load_orbitals()
// {
// psio_->open(PSIF_SCF_MOS,PSIO_OPEN_OLD);

// int basislength, old_puream;
// psio_->read_entry(PSIF_SCF_MOS,"BASIS NAME LENGTH",
// (char *)(&basislength),sizeof(int));
// char *basisnamec = new char[basislength];
// psio_->read_entry(PSIF_SCF_MOS,"BASIS NAME",basisnamec,
// basislength*sizeof(char));
// psio_->read_entry(PSIF_SCF_MOS,"PUREAM",(char *)(&old_puream),
// sizeof(int));
// bool old_forced_puream = (old_puream) ? true : false;
// std::string basisname(basisnamec);

// if (basisname == "")
// throw PSIEXCEPTION("SCF::load_orbitals: Custom basis sets not allowed for projection from a previous SCF");

// if (print_) {
// if (basisname != options_.get_str("BASIS")) {
// outfile->Printf(" Computing basis set projection from %s to %s.\n", \
// basisname.c_str(),options_.get_str("BASIS").c_str());
// } else {
// outfile->Printf(" Using orbitals from previous SCF, no projection.\n");
// }
// }

// std::shared_ptr<BasisSet> dual_basis;
// if (basisname != options_.get_str("BASIS")) {
// //std::shared_ptr<BasisSetParser> parser(new Gaussian94BasisSetParser(old_forced_puream));
// //molecule_->set_basis_all_atoms(basisname, "DUAL_BASIS_SCF");
// //dual_basis = BasisSet::construct(parser, molecule_, "DUAL_BASIS_SCF");
// dual_basis = BasisSet::pyconstruct_orbital(molecule_,
// "DUAL_BASIS_SCF", basisname, old_forced_puream);
// } else {
// dual_basis = BasisSet::pyconstruct_orbital(molecule_,
// "BASIS", options_.get_str("BASIS"));
// // TODO: oh my, forced_puream!
// // TODO: oh my, a basis for which a fn hasn't been set in the input translation
// // TODO: oh my, a non-fitting basis to be looked up (in Mol) not under BASIS
// //std::shared_ptr<BasisSet> dual_basis = BasisSet::pyconstruct(molecule_, basisname,
// // "DUAL_BASIS_SCF");
// // TODO: I think Rob was planning to rework this projection bit anyways
// // 2 Apr 2015: I (LAB) was hoping to avoid detangling this, but the need to
// // optimize w/custom basis sets has arrived before the new scf code, hence this hack
// }

// psio_->read_entry(PSIF_SCF_MOS,"SCF ENERGY",(char *) &(E_),sizeof(double));

// int old_nirrep, old_nsopi[8];
// psio_->read_entry(PSIF_SCF_MOS,"NIRREP",(char *) &(old_nirrep),sizeof(int));

// if (old_nirrep != nirrep_)
// throw PSIEXCEPTION("SCF::load_orbitals: Projection of orbitals between different symmetries is not currently supported");

// psio_->read_entry(PSIF_SCF_MOS,"NSOPI",(char *) (old_nsopi),nirrep_*sizeof(int));

// // Save current alpha and beta occupation vectors
// Dimension nalphapi_current (nirrep_, "Current number of alpha electrons per irrep");
// Dimension nbetapi_current (nirrep_, "Current number of beta electrons per irrep");
// nalphapi_current = nalphapi_;
// nbetapi_current = nbetapi_;

// // Read in guess alpha and beta occupation vectors
// psio_->read_entry(PSIF_SCF_MOS,"NALPHAPI",(char *) &(nalphapi_[0]),nirrep_*sizeof(int));
// psio_->read_entry(PSIF_SCF_MOS,"NBETAPI",(char *) &(nbetapi_[0]),nirrep_*sizeof(int));

// // Check if guess is the broken symmetry solution
// int nalpha_guess = 0;
// int nbeta_guess = 0;
// for (int h = 0; h < nirrep_; h++) {
// nalpha_guess += nalphapi_[h];
// nbeta_guess += nbetapi_[h];
// }

// bool guess_broken_symmetry = false;
// if (nalpha_guess == nbeta_guess && multiplicity_ == 3) guess_broken_symmetry = true;

// if (!broken_symmetry_) {
// if (!guess_broken_symmetry) {
// outfile->Printf( " Recomputing DOCC and SOCC from number of alpha and beta electrons from previous calculation.\n");
// for (int h = 0; h < nirrep_; h++) {
// doccpi_[h] = std::min(nalphapi_[h] , nbetapi_[h]);
// soccpi_[h] = std::abs(nalphapi_[h] - nbetapi_[h]);
// }
// print_occupation();

// SharedMatrix Ctemp_a(new Matrix("ALPHA MOS", nirrep_, old_nsopi, nalphapi_));
// Ctemp_a->load(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// SharedMatrix Ca;
// if (basisname != options_.get_str("BASIS")) {
// Ca = basis_projection(Ctemp_a, nalphapi_, dual_basis, basisset_);
// } else {
// Ca = Ctemp_a;
// }
// for (int h = 0; h < nirrep_; h++)
// for (int m = 0; m<nsopi_[h]; m++)
// for (int i = 0; i<nalphapi_[h]; i++)
// Ca_->set(h,m,i,Ca->get(h,m,i));

// SharedMatrix Ctemp_b(new Matrix("BETA MOS", nirrep_, old_nsopi, nbetapi_));
// Ctemp_b->load(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// SharedMatrix Cb;
// if (basisname != options_.get_str("BASIS")) {
// Cb = basis_projection(Ctemp_b, nbetapi_, dual_basis, basisset_);
// } else {
// Cb = Ctemp_b;
// }
// for (int h = 0; h < nirrep_; h++)
// for (int m = 0; m<nsopi_[h]; m++)
// for (int i = 0; i<nbetapi_[h]; i++)
// Cb_->set(h,m,i,Cb->get(h,m,i));
// }
// // If it's a triplet and there is a broken symmetry singlet guess - ignore the guess
// else{
// // Restore nalphapi and nbetapi requested by the user
// nalphapi_ = nalphapi_current;
// nbetapi_ = nbetapi_current;
// }
// }
// // if it's a broken symmetry solution
// else {

// SharedMatrix Ctemp_a(new Matrix("ALPHA MOS", nirrep_, old_nsopi, nalphapi_));
// Ctemp_a->load(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// SharedMatrix Ca;
// if (basisname != options_.get_str("BASIS")) {
// Ca = basis_projection(Ctemp_a, nalphapi_, dual_basis, basisset_);
// } else {
// Ca = Ctemp_a;
// }

// SharedMatrix Ctemp_b(new Matrix("BETA MOS", nirrep_, old_nsopi, nbetapi_));
// Ctemp_b->load(psio_, PSIF_SCF_MOS, Matrix::SubBlocks);
// SharedMatrix Cb;
// if (basisname != options_.get_str("BASIS")) {
// Cb = basis_projection(Ctemp_b, nbetapi_, dual_basis, basisset_);
// } else {
// Cb = Ctemp_b;
// }

// // Restore nalphapi and nbetapi requested by the user
// nalphapi_ = nalphapi_current;
// nbetapi_ = nbetapi_current;

// int socc_count = 0;
// for (int h = 0; h < nirrep_; h++) {
// // Copy doubly occupied orbitals into the orbital space for alpha and beta
// for (int i = 0; i<doccpi_[h]; i++) {
// for (int m = 0; m<nsopi_[h]; m++) {
// Ca_->set(h,m,i,Ca->get(h,m,i));
// Cb_->set(h,m,i,Cb->get(h,m,i));
// }
// }
// // Copy singly occupied orbitals into the appropriate alpha and beta orbital spaces
// for (int i = doccpi_[h]; i<(doccpi_[h]+soccpi_[h]); i++) {
// socc_count++;
// if (socc_count == 1) {
// for (int m = 0; m<nsopi_[h]; m++) {
// Ca_->set(h,m,doccpi_[h],Ca->get(h,m,i));
// }
// }
// if (socc_count == 2) {
// for (int m = 0; m<nsopi_[h]; m++) {
// Cb_->set(h,m,doccpi_[h],Ca->get(h,m,i));
// }
// }
// }
// }
// }
// psio_->close(PSIF_SCF_MOS,1);
// delete[] basisnamec;
// }
void HF::format_guess()
{
// Nothing to do, only for special cases
}


void HF::check_phases()
Expand Down Expand Up @@ -2091,6 +1829,8 @@ void HF::print_energies()
Process::environment.globals["DISPERSION CORRECTION ENERGY"] = energies_["-D"];
}

Process::environment.globals["SCF ITERATIONS"] = iteration_;

// Only print this alert if we are actually doing EFP or PCM
if(pcm_enabled_ || ( Process::environment.get_efp()->get_frag_count() > 0 ) ) {
outfile->Printf(" Alert: EFP and PCM quantities not currently incorporated into SCF psivars.");
Expand Down
3 changes: 3 additions & 0 deletions psi4/src/psi4/libscf_solver/hf.h
View file
Open in desktop
Expand Up @@ -393,6 +393,9 @@ class HF : public Wavefunction {
/** Form X'(FDS - SDF)X (for DIIS) **/
virtual SharedMatrix form_FDSmSDF(SharedMatrix Fso, SharedMatrix Dso);

/** Performs any operations required for a incoming guess **/
virtual void format_guess();

/** Save orbitals to use later as a guess **/
// virtual void save_orbitals();

Expand Down
17 changes: 16 additions & 1 deletion psi4/src/psi4/libscf_solver/rohf.cc
View file
Open in desktop
Expand Up @@ -78,7 +78,7 @@ void ROHF::common_init()
moFeff_ = SharedMatrix(factory_->create_matrix("F effective (MO basis)"));
soFeff_ = SharedMatrix(factory_->create_matrix("F effective (orthogonalized SO basis)"));
Ct_ = SharedMatrix(factory_->create_matrix("Orthogonalized Molecular orbitals"));
Ca_ = SharedMatrix(factory_->create_matrix("C"));
Ca_ = SharedMatrix(factory_->create_matrix("C"));
Cb_ = Ca_;
Da_ = SharedMatrix(factory_->create_matrix("SCF alpha density"));
Db_ = SharedMatrix(factory_->create_matrix("SCF beta density"));
Expand All @@ -104,6 +104,20 @@ void ROHF::common_init()
}
}

void ROHF::format_guess()
{
// Need to build Ct_

// Canonical Orthogonalization
if (X_->rowspi() != X_->colspi()){
throw PSIEXCEPTION("ROHF::format_guess: 'GUESS READ' is not available for canonical orthogonalization cases.");

// Symmetric Orthogonalization
} else {
Ct_ = Matrix::triplet(X_, S_, Ca_);
}
}

void ROHF::semicanonicalize()
{
// Quick sanity checks
Expand Down Expand Up @@ -376,6 +390,7 @@ void ROHF::form_initialF()

void ROHF::form_F()
{

// Start by constructing the standard Fa and Fb matrices encountered in UHF
Fa_->copy(H_);
Fb_->copy(H_);
Expand Down
1 change: 1 addition & 0 deletions psi4/src/psi4/libscf_solver/rohf.h
View file
Open in desktop
Expand Up @@ -73,6 +73,7 @@ class ROHF : public HF {
virtual void finalize();

void save_density_and_energy();
void format_guess();

// Second-order convergence code
void Hx(SharedMatrix x, SharedMatrix ret);
Expand Down
2 changes: 1 addition & 1 deletion tests/CMakeLists.txt
View file
Open in desktop
Expand Up @@ -78,7 +78,7 @@ foreach(test_name adc1 adc2 casscf-fzc-sp casscf-sa-sp ao-casscf-sp casscf-sp ca
pywrap-db3 pywrap-freq-e-sowreap pywrap-freq-g-sowreap
pywrap-molecule pywrap-opt-sowreap rasci-c2-active rasci-h2o
rasci-ne rasscf-sp sad1 sapt1 sapt2 sapt3 sapt4 sapt5 sapt6
sapt7 sapt8 scf-bz2 scf-guess-read scf-bs scf1 scf-occ
sapt7 sapt8 scf-bz2 scf-guess-read1 scf-guess-read2 scf-bs scf1 scf-occ
scf2 scf3 scf4 scf5 scf6 scf-property soscf1 soscf2 stability1
stability2 tu1-h2o-energy tu2-ch2-energy tu3-h2o-opt
tu4-h2o-freq tu5-sapt tu6-cp-ne2 x2c1 x2c2 x2c3 zaptn-nh2
Expand Down
3 changes: 0 additions & 3 deletions tests/scf-guess-read/CMakeLists.txt
View file
Open in desktop

This file was deleted.

3 changes: 3 additions & 0 deletions tests/scf-guess-read1/CMakeLists.txt
View file
Open in desktop
@@ -0,0 +1,3 @@
include(TestingMacros)

add_regression_test(scf-guess-read1 "psi;quicktests;scf")
0 tests/scf-guess-read/input.dat → tests/scf-guess-read1/input.dat
View file
Open in desktop
File renamed without changes.
0 tests/scf-guess-read/output.ref → tests/scf-guess-read1/output.ref
View file
Open in desktop
File renamed without changes.
3 changes: 3 additions & 0 deletions tests/scf-guess-read2/CMakeLists.txt
View file
Open in desktop
@@ -0,0 +1,3 @@
include(TestingMacros)

add_regression_test(scf-guess-read2 "psi;scf")