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

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

Merged
merged 11 commits into from Jul 3, 2017

Conversation

Projects
None yet
4 participants
@ashutoshvt
Member

ashutoshvt commented Jun 26, 2017

Description

This pull request aims to properly update energy, orbital energies, densities etc of the wave function for (conventional and density-fitted Mp2) and density-fitted CCSD calculations, so as to ensure proper checkpointing using FCHK interface.
Also, there seems to be a problem in using FCHKWriter for conventional mp2 gradient calculations as I get a segmentation fault. This PR is meant to fix these problems which was pointed out by abdale on the psi4 forum. http://forum.psicode.org/t/problems-with-wavefunction-object-and-fchk-interface/532

Status

  • Ready to go
@dgasmith

Thanks! The rest of it looks fr at however, I am starting to worry about how many post-SCF methods actually set their 'energy_'.

Ca_ = SharedMatrix(new Matrix("DF-MP2 (unrelaxed) Natural Orbitals", nsopi_, nmopi_));
epsilon_a_ = SharedVector(new Vector("DF-MP2 (unrelaxed) NO Occupations", nmopi_));
Ca__ = SharedMatrix(new Matrix("DF-MP2 (unrelaxed) Natural Orbitals", nsopi_, nmopi_));
epsilon_a__ = SharedVector(new Vector("DF-MP2 (unrelaxed) NO Occupations", nmopi_));

This comment has been minimized.

@dgasmith

dgasmith Jun 27, 2017

Member

I would think that we actually want the non relaxed natural orbitals for better guesses downstream. I wouldn't add a Ca__ attribute as it's difficult to tell between Ca_ and Ca__.

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 27, 2017

Member

As you pointed out, many of the post-scf methods don't set their energy_ and just update the process environment variable. This needs a major overhaul. Here, I just wanted to fix the problem faced by a user on the forum. Also, I was uncertain if we needed the unrelaxed natural orbitals coefficients and occupation numbers as I didn't see them getting used anywhere. How do you suggest I retain these as I want the Ca_ and epsilon_a_ to be the conventional MO coefficients and orbital energies respectively.

This comment has been minimized.

@dgasmith

dgasmith Jun 27, 2017

Member

I don't think we want to retain Ca_ and epsilon_a_ as the original SCF values, its ok if Wavefunctions update themselves. For example, the MCSCF Wavefunction rotates its own orbitals.

These values will not get used within the MP2 code, but are useful for things like:

set mp2_opdm true # or whatever
mp2_grad, mp2_wfn = psi4.energy("MP2")
psi4.energy("OCCSD", ref_wfn = mp2_wfn)

The relaxed densities are odd because they rotate the occupied x virtual parts of the OPDM while the unrelaxed density only rotates the occupied x occupied and virtual x virtual areas. The former actually changes the Wavefunction while the latter helps compress the information. (Hopefully I didn't screw this up, its been awhile since I have done MP2 densities)

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 27, 2017

Member

I see. However, its confusing when one uses the fchk_interface as Ca_ and epsilon_a_ which are natural orbitals and occupation numbers get printed as Alpha MO coefficients and Alpha Orbital Energies respectively in the checkpoint file. Do you think we should provide the original MO coefficients and orbital energies of the reference wavefunction as well in the df-mp2 gradient case since all the other modules do the same?

This comment has been minimized.

@dgasmith

dgasmith Jun 27, 2017

Member

I don't think all other modules provide the original MO coefficients, in fact in many cases they alter the coefficients as described above. If someone wants SCF level fchk, they should give the fchk a SCF Wavefunction? We want the fchk written to write the current state of the Wavefunction, regardless of what that state is I would think.

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 27, 2017

Member

You are right, I was thinking of including the SCF Wavefunction data for post-scf methods which i guess can be done quite easily from the python's side:

scf_e, scf_wfn = energy('scf', return_wfn=True)
scf_fchk_writer = core.FCHKWriter(scf_wfn)
scf_fchk_writer.write('scf.fchk')
mp2_e, mp2_wfn = gradient('mp2', ref_wfn = scf_wfn, return_wfn=True)
mp2_fchk_writer = core.FCHKWriter(mp2_wfn)
mp2_fchk_writer.write('mp2.fchk')

Also, do you think we should change the name Alpha Orbital Energies which actually are the occupation numbers in the case of df-mp2 to something more suitable?

This comment has been minimized.

@dgasmith

dgasmith Jun 27, 2017

Member

Good question. @andysim @loriab I forget who wrote the fchk interface. I unfortunately do not know a lot about the spec.

This comment has been minimized.

@loriab

loriab Jun 27, 2017

Member

Its not unheard of to have occ num in place of energies. That's how I confirmed was dealing with natorbs in cfour. There does seem to be an Occup slot in Molden spec, though.

 Sym= symmetry_label_1
 Ene= mo_energy_1
 Spin= (Alpha|Beta)
 Occup= mo_occupation_number_1
 ao_number_1 mo_coefficient_1
 ...
 ao_number_n mo_coefficient_n
 ....

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 27, 2017

Member

Thanks for the reply @loriab I was wondering if we can print the name of the matrix on the fchk file, something like
wfn_->Ca()->name()
instead of hardcoding it like "Alpha Orbital Energies"
write_matrix("Alpha Orbital Energies", wavefunction_->epsilon_a_subset("AO"));
This way it will reflect whatever change was made by the corresponding wavefunction ?

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 27, 2017

Member

I will push appropriate changes in this regard which would make it easier to get feedback from you guys.

@ashutoshvt

This comment has been minimized.

Member

ashutoshvt commented Jun 28, 2017

@dgasmith I think this is ready for feedback.

@@ -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)"));

This comment has been minimized.

@dgasmith

dgasmith Jun 28, 2017

Member

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

This comment has been minimized.

@ashutoshvt

ashutoshvt Jun 28, 2017

Member

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?

@loriab

loriab approved these changes Jun 30, 2017

lgtm

@@ -33,6 +33,8 @@
#include "defines.h"
#include "psi4/libpsi4util/PsiOutStream.h"
#include "psi4/liboptions/liboptions.h"
#include "psi4/libpsi4util/process.h"

This comment has been minimized.

@loriab

loriab Jun 30, 2017

Member

Should this header go in dfocc/manager.cc so as to minimize its impact on the build?

@ashutoshvt

This comment has been minimized.

Member

ashutoshvt commented Jun 30, 2017

@loriab Yup, removed that header now.

@dgasmith

This comment has been minimized.

Member

dgasmith commented Jul 1, 2017

LGTM @andysim Can you check this one out, I think you may know a bit about fchk.

@@ -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;

This comment has been minimized.

@andysim

andysim Jul 3, 2017

Member

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

This comment has been minimized.

@ashutoshvt

ashutoshvt Jul 3, 2017

Member

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

@andysim

andysim approved these changes Jul 3, 2017

Looks great, thanks! I hacked the original FCHK stuff together, based entirely on the few files that ship with GDMA, so I'm no expert at all. I think it's possible to have both the correlated and SCF densities in there, but don't have any way to test that (I can't see a standard for the file format in my quick Google search). Regardless of whether we have that feature, this is a much needed patch.

@ashutoshvt

This comment has been minimized.

Member

ashutoshvt commented Jul 3, 2017

I agree, there are no well-defined standards as of yet. @dgasmith suggested that one can do a scf calculation to get scf level properties as this info is meant for the state of the wfn requested.

@dgasmith

This comment has been minimized.

Member

dgasmith commented Jul 3, 2017

I think in terms of changes for this particular patch this is all ready to go. Any objections?

@ashutoshvt

This comment has been minimized.

Member

ashutoshvt commented Jul 3, 2017

No objections from my side.

@andysim andysim merged commit 00d0f83 into psi4:master Jul 3, 2017

2 checks passed

continuous-integration/Distelli
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@loriab loriab added this to the Psi4 1.2 milestone Apr 13, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment