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.

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

Proper update of post-scf wave-functions for fchk interface #747

Merged
merged 11 commits into from Jul 3, 2017
Merged

Proper update of post-scf wave-functions for fchk interface #747

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
a725aa8
updated energy of df-mp2 wavefunction
ashutoshvt Jun 8, 2017
d63e619
updated Mp2 energy for occ and dfocc wfns
ashutoshvt Jun 25, 2017
8c53bd6
prevents change in Ca_ and epsilon_a in dfmp2
ashutoshvt Jun 25, 2017
f1b5dd6
cleaned up and updated names of wfns
ashutoshvt Jun 25, 2017
e416afe
modified names of orbital energies and coefficients in fchk
ashutoshvt Jun 27, 2017
b84fad4
Ca_ is not reset now for future use
ashutoshvt Jun 27, 2017
e3b66d6
modified names of coefficients and orbital energies in hf,
ashutoshvt Jun 27, 2017
a0f1ce9
updated energies and names for ccsd wfns
ashutoshvt Jun 28, 2017
f513e93
removed unnecessary header
ashutoshvt Jun 30, 2017
1beda8e
updated name and energy ccsd(t) wfn
ashutoshvt Jul 3, 2017
3dc70e9
updated names and energies of all wfns in dfocc/manager.cc
ashutoshvt Jul 3, 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
1 change: 1 addition & 0 deletions psi4/src/psi4/ccenergy/ccenergy.cc
View file
Open in desktop
Expand Up @@ -514,6 +514,7 @@ double CCEnergyWavefunction::compute_energy()
#endif

energy_ = moinfo_.ecc + moinfo_.eref;
name_ = "CCSD";
Process::environment.globals["CURRENT ENERGY"] = moinfo_.ecc+moinfo_.eref;
Process::environment.globals["CURRENT CORRELATION ENERGY"] = moinfo_.ecc;
//Process::environment.globals["CC TOTAL ENERGY"] = moinfo.ecc+moinfo.eref;
Expand Down
7 changes: 7 additions & 0 deletions psi4/src/psi4/dfmp2/mp2.cc
View file
Open in desktop
Expand Up @@ -229,6 +229,7 @@ double DFMP2::compute_energy()
form_energy();
timer_off("DFMP2 Energy");
print_energies();
energy_ = variables_["MP2 TOTAL ENERGY"];

return variables_["MP2 TOTAL ENERGY"];
}
Expand Down Expand Up @@ -314,6 +315,7 @@ SharedMatrix DFMP2::compute_gradient()

if(options_.get_bool("ONEPDM")){
print_energies();
energy_ = variables_["MP2 TOTAL ENERGY"];
return SharedMatrix(new Matrix("NULL", 0, 0));
}

Expand All @@ -322,6 +324,7 @@ SharedMatrix DFMP2::compute_gradient()
timer_off("DFMP2 grad");

print_energies();
energy_ = variables_["MP2 TOTAL ENERGY"];

print_gradients();

Expand Down Expand Up @@ -2307,6 +2310,7 @@ void RDFMP2::form_Z()
SharedMatrix AP(new Matrix("A_mn^ls P_ls^(2)", nso, nso));
double** APp = AP->pointer();
SharedMatrix Dtemp;


if(options_.get_bool("OPDM_RELAX")){
psio_->read_entry(PSIF_DFMP2_AIA, "W", (char*) Wpq1p[0], sizeof(double) * nmo * nmo);
Expand Down Expand Up @@ -2403,6 +2407,7 @@ void RDFMP2::form_Z()
Ca_ = SharedMatrix(new Matrix("DF-MP2 Natural Orbitals", nsopi_, nmopi_));
epsilon_a_ = SharedVector(new Vector("DF-MP2 NO Occupations", nmopi_));
Da_ = SharedMatrix(new Matrix("DF-MP2 relaxed density", nsopi_, nsopi_));

}else{
// Don't relax the OPDM
Dtemp = Ppq->clone();
Expand All @@ -2411,11 +2416,13 @@ void RDFMP2::form_Z()
Dtemp->scale(0.5);

// Add in the reference contribution

for (int i = 0; i < nocc; ++i)
Dtemp->add(i, i, 1.0);
Ca_ = SharedMatrix(new Matrix("DF-MP2 (unrelaxed) Natural Orbitals", nsopi_, nmopi_));
epsilon_a_ = SharedVector(new Vector("DF-MP2 (unrelaxed) NO Occupations", nmopi_));
Da_ = SharedMatrix(new Matrix("DF-MP2 unrelaxed density", nsopi_, nsopi_));

}

compute_opdm_and_nos(Dtemp, Da_, Ca_, epsilon_a_);
Expand Down
1 change: 1 addition & 0 deletions psi4/src/psi4/dfocc/dfocc.cc
View file
Open in desktop
Expand Up @@ -35,6 +35,7 @@
#include "psi4/liboptions/liboptions.h"



using namespace psi;


Expand Down
36 changes: 36 additions & 0 deletions psi4/src/psi4/dfocc/manager.cc
View file
Open in desktop
Expand Up @@ -511,6 +511,11 @@ void DFOCC::mp2_manager()
Process::environment.globals["MP2 OPPOSITE-SPIN CORRELATION ENERGY"] = Emp2AB;
Process::environment.globals["MP2 SAME-SPIN CORRELATION ENERGY"] = Emp2AA+Emp2BB;

/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["MP2 TOTAL ENERGY"];
name_ = "DF-MP2";


// S2
//if (comput_s2_ == "TRUE" && reference_ == "UNRESTRICTED" && dertype == "NONE") s2_response();
if (comput_s2_ == "TRUE" && reference_ == "UNRESTRICTED") {
Expand Down Expand Up @@ -804,6 +809,10 @@ void DFOCC::ccsd_manager()
Process::environment.globals["CCSD TOTAL ENERGY"] = Eccsd;
Process::environment.globals["CCSD CORRELATION ENERGY"] = Eccsd - Escf;

/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCSD TOTAL ENERGY"];
name_ = "DF-CCSD";

// CCSDL
if (dertype == "FIRST" || cc_lambda_ == "TRUE") {
// memalloc
Expand Down Expand Up @@ -1171,6 +1180,9 @@ void DFOCC::ccsd_t_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = Eccsd_t - Escf;
Process::environment.globals["CCSD(T) TOTAL ENERGY"] = Eccsd_t;
Process::environment.globals["(T) CORRECTION ENERGY"] = E_t;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCSD(T) TOTAL ENERGY"];
name_ = "DF-CCSD(T)";

/*
// CCSDL
Expand Down Expand Up @@ -1551,6 +1563,9 @@ void DFOCC::ccsdl_t_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = Eccsd_at - Escf;
Process::environment.globals["CCSD(AT) TOTAL ENERGY"] = Eccsd_at;
Process::environment.globals["(AT) CORRECTION ENERGY"] = E_at;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCSD(AT) TOTAL ENERGY"];
name_ = "DF-CCSD(AT)";

/*
// Compute Analytic Gradients
Expand Down Expand Up @@ -1816,6 +1831,9 @@ void DFOCC::ccd_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = Eccd - Escf;
Process::environment.globals["CCD TOTAL ENERGY"] = Eccd;
Process::environment.globals["CCD CORRELATION ENERGY"] = Eccd - Escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCD TOTAL ENERGY"];
name_ = "DF-CCD";

// CCDL
if (dertype == "FIRST" || cc_lambda_ == "TRUE") {
Expand Down Expand Up @@ -2141,6 +2159,9 @@ void DFOCC::omp3_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = Emp3L - Escf;
Process::environment.globals["OMP3 TOTAL ENERGY"] = Emp3L;
Process::environment.globals["OMP3 CORRELATION ENERGY"] = Emp3L - Escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["OMP3 TOTAL ENERGY"];
name_ = "DF-OMP3";

// OEPROP
if (oeprop_ == "TRUE") oeprop();
Expand Down Expand Up @@ -2356,6 +2377,9 @@ void DFOCC::mp3_manager()
Process::environment.globals["MP3 TOTAL ENERGY"] = Emp3;
Process::environment.globals["MP3 CORRELATION ENERGY"] = Emp3 - Escf;
Emp3L=Emp3;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["MP3 TOTAL ENERGY"];
name_ = "DF-MP3";

// Compute Analytic Gradients
if (dertype == "FIRST" || ekt_ip_ == "TRUE") {
Expand Down Expand Up @@ -2665,6 +2689,9 @@ void DFOCC::omp2_5_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = Emp3L - Escf;
Process::environment.globals["OMP2.5 TOTAL ENERGY"] = Emp3L;
Process::environment.globals["OMP2.5 CORRELATION ENERGY"] = Emp3L - Escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["OMP2.5 TOTAL ENERGY"];
name_ = "DF-OMP2.5";

// OEPROP
if (oeprop_ == "TRUE") oeprop();
Expand Down Expand Up @@ -2861,6 +2888,9 @@ void DFOCC::mp2_5_manager()
Process::environment.globals["MP2.5 TOTAL ENERGY"] = Emp3;
Process::environment.globals["MP2.5 CORRELATION ENERGY"] = Emp3 - Escf;
Emp3L=Emp3;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["MP2.5 TOTAL ENERGY"];
name_ = "DF-MP2.5";

// Compute Analytic Gradients
if (dertype == "FIRST" || ekt_ip_ == "TRUE") {
Expand Down Expand Up @@ -3143,6 +3173,9 @@ void DFOCC::olccd_manager()
Process::environment.globals["CURRENT CORRELATION ENERGY"] = ElccdL - Escf;
Process::environment.globals["OLCCD TOTAL ENERGY"] = ElccdL;
Process::environment.globals["OLCCD CORRELATION ENERGY"] = ElccdL - Escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["OLCCD TOTAL ENERGY"];
name_ = "DF-OLCCD";

// OEPROP
if (oeprop_ == "TRUE") oeprop();
Expand Down Expand Up @@ -3336,6 +3369,9 @@ void DFOCC::lccd_manager()
Process::environment.globals["LCCD TOTAL ENERGY"] = Elccd;
Process::environment.globals["LCCD CORRELATION ENERGY"] = Elccd - Escf;
ElccdL = Elccd;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["LCCD TOTAL ENERGY"];
name_ = "DF-LCCD";

// Compute Analytic Gradients
if (dertype == "FIRST" || ekt_ip_ == "TRUE") {
Expand Down
3 changes: 3 additions & 0 deletions psi4/src/psi4/fnocc/ccsd.cc
View file
Open in desktop
Expand Up @@ -201,6 +201,9 @@ double CoupledCluster::compute_energy() {
Process::environment.globals["CCSD OPPOSITE-SPIN CORRELATION ENERGY"] = eccsd_os;
Process::environment.globals["CCSD SAME-SPIN CORRELATION ENERGY"] = eccsd_ss;
Process::environment.globals["CCSD TOTAL ENERGY"] = eccsd + escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCSD TOTAL ENERGY"];
name_ = "CCSD";
}else{
Process::environment.globals["QCISD CORRELATION ENERGY"] = eccsd;
Process::environment.globals["QCISD OPPOSITE-SPIN CORRELATION ENERGY"] = eccsd_os;
Expand Down
3 changes: 3 additions & 0 deletions psi4/src/psi4/fnocc/df_ccsd.cc
View file
Open in desktop
Expand Up @@ -102,6 +102,9 @@ double DFCoupledCluster::compute_energy() {
Process::environment.globals["CCSD SAME-SPIN CORRELATION ENERGY"] = eccsd_ss;
Process::environment.globals["CCSD TOTAL ENERGY"] = eccsd + escf;
Process::environment.globals["CURRENT ENERGY"] = eccsd + escf;
/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["CCSD TOTAL ENERGY"];
name_ = "DF-CCSD";

if (options_.get_bool("COMPUTE_TRIPLES")){
long int o = ndoccact;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we get in here, shouldn't the name be updated to DF-CCSD(T)?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks for pointing that out @andysim. I took this opportunity to update the names and energies of all the wfns inside manager.cc

Expand Down
5 changes: 4 additions & 1 deletion psi4/src/psi4/fnocc/fnocc.cc
View file
Open in desktop
Expand Up @@ -72,9 +72,11 @@ SharedWavefunction fnocc(SharedWavefunction ref_wfn, Options &options) {
if ( !options.get_bool("RUN_CEPA") ) {
std::shared_ptr<CoupledCluster> ccsd(new CoupledCluster(wfn,options));
ccsd->compute_energy();
return ccsd;
} else {
std::shared_ptr<CoupledPair> cepa (new CoupledPair(wfn,options));
cepa->compute_energy();
return cepa;
}

}else {
Expand Down Expand Up @@ -114,10 +116,11 @@ SharedWavefunction fnocc(SharedWavefunction ref_wfn, Options &options) {

tstop();

return ccsd;

}

return wfn;
//return wfn;
} // end fnocc

}} // end namespaces
3 changes: 3 additions & 0 deletions psi4/src/psi4/libmints/matrix.cc
View file
Open in desktop
Expand Up @@ -91,6 +91,7 @@ Matrix::Matrix(const Matrix &c)
matrix_ = NULL;
nirrep_ = c.nirrep_;
symmetry_ = c.symmetry_;
name_ = c.name();
alloc();
copy_from(c.matrix_);
}
Expand All @@ -101,6 +102,7 @@ Matrix::Matrix(const SharedMatrix &c)
matrix_ = NULL;
nirrep_ = c->nirrep_;
symmetry_ = c->symmetry_;
name_ = c->name();
alloc();
copy_from(c->matrix_);
}
Expand All @@ -111,6 +113,7 @@ Matrix::Matrix(const Matrix *c)
matrix_ = NULL;
nirrep_ = c->nirrep_;
symmetry_ = c->symmetry_;
name_ = c->name();
alloc();
copy_from(c->matrix_);
}
Expand Down
15 changes: 10 additions & 5 deletions psi4/src/psi4/libmints/writer.cc
View file
Open in desktop
Expand Up @@ -78,7 +78,6 @@ void GradientWriter::write(const std::string &filename)
MoldenWriter::MoldenWriter(std::shared_ptr<Wavefunction> wavefunction)
: wavefunction_(wavefunction)
{

}
void MoldenWriter::write(const std::string &filename, std::shared_ptr<Matrix> Ca, std::shared_ptr<Matrix> Cb, std::shared_ptr<Vector> Ea, std::shared_ptr<Vector> Eb, std::shared_ptr<Vector> OccA, std::shared_ptr<Vector> OccB, bool dovirtual)
{
Expand Down Expand Up @@ -274,6 +273,8 @@ void MoldenWriter::write(const std::string &filename, std::shared_ptr<Matrix> Ca
FCHKWriter::FCHKWriter(std::shared_ptr<Wavefunction> wavefunction)
: wavefunction_(wavefunction)
{
SharedMatrix Ca = wavefunction_->Ca();
Ca->print();
}


Expand Down Expand Up @@ -634,10 +635,14 @@ void FCHKWriter::write(const std::string &filename)
write_matrix("Contraction coefficients", coefficients);
write_matrix("Coordinates of each shell", shell_coords);
write_number("Total Energy", wavefunction_->reference_energy());
write_matrix("Alpha Orbital Energies", wavefunction_->epsilon_a_subset("AO"));
write_matrix("Alpha MO coefficients", reorderedCa);
write_matrix("Beta Orbital Energies", wavefunction_->epsilon_b_subset("AO"));
write_matrix("Beta MO coefficients", reorderedCb);
//write_matrix("Alpha Orbital Energies", wavefunction_->epsilon_a_subset("AO"));
write_matrix(wavefunction_->epsilon_a()->name().c_str(), wavefunction_->epsilon_a_subset("AO"));
//write_matrix("Alpha MO coefficients", reorderedCa);
write_matrix(wavefunction_->Ca()->name().c_str(), reorderedCa);
//write_matrix("Beta Orbital Energies", wavefunction_->epsilon_b_subset("AO"));
write_matrix(wavefunction_->epsilon_b()->name().c_str(), wavefunction_->epsilon_b_subset("AO"));
//write_matrix("Beta MO coefficients", reorderedCb);
write_matrix(wavefunction_->Cb()->name().c_str(), reorderedCb);
char* label = new char[256];
std::string type = name == "DFT" ? "SCF" : name;
sprintf(label, "Total %s Density", type.c_str());
Expand Down
3 changes: 2 additions & 1 deletion psi4/src/psi4/libscf_solver/rhf.cc
View file
Open in desktop
Expand Up @@ -87,9 +87,10 @@ void RHF::common_init()
// Allocate matrix memory
Fa_ = SharedMatrix(factory_->create_matrix("F"));
Fb_ = Fa_;
Ca_ = SharedMatrix(factory_->create_matrix("C"));
Ca_ = SharedMatrix(factory_->create_matrix("MO coefficients (C)"));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these new names must be identical to that of the replaced names in writer.cc for this for fichk to work?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fchk will work even if the names were different. For example, If one looks at mp2.cc inside dfmp2, Ca_ = SharedMatrix(new Matrix("DF-MP2 Natural Orbitals", nsopi_, nmopi_));
I wanted this to be reflected when we write the wfn's info to the fchk_writer If the wfn doesn't do anything special, the names would stay the same as reference wfn's. Does this make sense?

Cb_ = Ca_;
epsilon_a_ = SharedVector(factory_->create_vector());
epsilon_a_->set_name("orbital energies");
epsilon_b_ = epsilon_a_;
Da_ = SharedMatrix(factory_->create_matrix("SCF density"));
Db_ = Da_;
Expand Down
5 changes: 4 additions & 1 deletion psi4/src/psi4/libscf_solver/rohf.cc
View file
Open in desktop
Expand Up @@ -83,8 +83,9 @@ 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("alpha MO coefficients (C)"));
Cb_ = Ca_;
Cb_->set_name("beta MO coefficients (C)");
Da_ = SharedMatrix(factory_->create_matrix("SCF alpha density"));
Db_ = SharedMatrix(factory_->create_matrix("SCF beta density"));
Lagrangian_ = SharedMatrix(factory_->create_matrix("Lagrangian matrix"));
Expand All @@ -100,7 +101,9 @@ void ROHF::common_init()
moFb_ = SharedMatrix(factory_->create_matrix("MO beta Fock Matrix (MO basis)"));

epsilon_a_ = SharedVector(factory_->create_vector());
epsilon_a_->set_name("alpha orbital energies");
epsilon_b_ = epsilon_a_;
epsilon_b_->set_name("beta orbital energies");
same_a_b_dens_ = false;
same_a_b_orbs_ = true;

Expand Down
6 changes: 4 additions & 2 deletions psi4/src/psi4/libscf_solver/uhf.cc
View file
Open in desktop
Expand Up @@ -95,8 +95,8 @@ void UHF::common_init()
Db_old_ = SharedMatrix(factory_->create_matrix("Old beta SCF density"));
Dt_old_ = SharedMatrix(factory_->create_matrix("D total old"));
Lagrangian_ = SharedMatrix(factory_->create_matrix("Lagrangian"));
Ca_ = SharedMatrix(factory_->create_matrix("C alpha"));
Cb_ = SharedMatrix(factory_->create_matrix("C beta"));
Ca_ = SharedMatrix(factory_->create_matrix("alpha MO coefficients (C)"));
Cb_ = SharedMatrix(factory_->create_matrix("beta MO coefficients (C)"));
Ga_ = SharedMatrix(factory_->create_matrix("G alpha"));
Gb_ = SharedMatrix(factory_->create_matrix("G beta"));
Va_ = SharedMatrix(factory_->create_matrix("V alpha"));
Expand All @@ -108,7 +108,9 @@ void UHF::common_init()
wKb_ = SharedMatrix(factory_->create_matrix("wK beta"));

epsilon_a_ = SharedVector(factory_->create_vector());
epsilon_a_->set_name("alpha orbital energies");
epsilon_b_ = SharedVector(factory_->create_vector());
epsilon_b_->set_name("beta orbital energies");

same_a_b_dens_ = false;
same_a_b_orbs_ = false;
Expand Down
5 changes: 5 additions & 0 deletions psi4/src/psi4/occ/manager.cc
View file
Open in desktop
Expand Up @@ -385,6 +385,11 @@ void OCCWave::mp2_manager()
}
}

/* updates the wavefunction for checkpointing */
energy_ = Process::environment.globals["MP2 TOTAL ENERGY"];
name_ = "MP2";


// S2
//if (comput_s2_ == "TRUE" && reference_ == "UNRESTRICTED") s2_response();

Expand Down
2 changes: 1 addition & 1 deletion psi4/src/psi4/occ/occwave.cc
View file
Open in desktop
Expand Up @@ -593,7 +593,7 @@ void OCCWave::mem_release()
delete oo_pairidxAA;
delete vv_pairidxAA;

Ca_.reset();
//Ca_.reset();
Ca_ref.reset();
Hso.reset();
Tso.reset();
Expand Down