DLPNO-MP2 #2093
DLPNO-MP2 #2093
Conversation
There was a problem hiding this comment.
yay for new functionality!
Your mp2_type/qc_module question is an incisive one. Generally, a new approach will show up in either the method name (energy(methodname-mp2)) or in the algorithm option (mp2_type=dlpno), not both. So algorithm options (disk_df vs. mem_df) are definitely mp2_type. broad categories that give different answers (conv, df, cd) have usually been formuated as mp2_type. methods that give different answers and are known by name (e.g., fno) have been forulated by method name. I can see dlpno going either way. But the most straightforward for now would be to continue following fno, and use energy("dlpno-mp2") and not define a qc_module or mp2_type for dlpno at present.
For conv options, it'd be good to tie those options (if not set explicitly) into the main e_conv, d/r_conv options.
In fnocc, AED allows setting some conv crit by percentage, energy, and cutoff. relevant here?
|
Fantastic! We want to hook this up to double-hybrid DFT eventually, then a Users familiar with ORCA will expect the Loose/Normal/TightPNO quick settings. |
Both good ideas! I added a |
Good to know the difference between module and type, I'll leave as Per Holger's suggestion, I've added a single DLPNO's |
psi4/src/psi4/dlpno/mp2.cc
Outdated
| // atomic mulliken populations for this orbital | ||
| std::vector<double> mkn_pop(natom, 0.0); | ||
|
|
||
| SharedMatrix P_i = reference_wavefunction_->S()->clone(); | ||
|
|
||
| for(size_t u=0; u < nbf; u++){ | ||
| P_i->scale_row(0, u, C_lmo->get(u, i)); | ||
| P_i->scale_column(0, u, C_lmo->get(u, i)); | ||
| } | ||
|
|
||
| for(size_t u=0; u < nbf; u++){ | ||
| int centerU = basisset_->function_to_center(u); | ||
| double p_uu = P_i->get(u,u); | ||
|
|
||
| for(size_t v=0; v < nbf; v++){ | ||
| int centerV = basisset_->function_to_center(u); | ||
| double p_vv = P_i->get(v,v); | ||
|
|
||
| // off-diag pops (p_uv) split between u and v prop to diag pops | ||
| double p_uv = P_i->get(u,v); | ||
| mkn_pop[centerU] += p_uv * ((p_uu) / (p_uu + p_vv)); | ||
| mkn_pop[centerV] += p_uv * ((p_vv) / (p_uu + p_vv)); | ||
|
|
||
| } | ||
| } |
There was a problem hiding this comment.
I'm wondering whether anyone has done this with better methods, e.g. IAO. Mulliken doesn't even have a basis set limit...
There was a problem hiding this comment.
The other part of the algorithms should work just the same, regardless of how the partial charge is computed.
There was a problem hiding this comment.
Turbomole's OSV-PNO uses IAOs. Maybe Molpro's variant as well.
ORCA had initially some problems with finding different LMO solutions, e.g. for augmented basis sets, and worked on getting a really robust localizer (augmented hessian algorithm).
So looking into what LMOs PSI4 should use is something worth spending time on. Long term, not this PR.
There was a problem hiding this comment.
IAOs are a simple solution; they're almost a drop-in replacement to Mulliken. One can also come up with other mathematically well-defined partial charges, as we discuss in https://pubs.acs.org/doi/10.1021/ct401016x
There was a problem hiding this comment.
That's a cool paper Susi, thanks. IAO charges would be an improvement over Mulliken charges; I know Molpro's related PNO-LMP2 method uses IAO charges for domain construction. There's some IAO/IBO code buried in the fisapt module. Moving that code to the globally accessible Localizer class would be a worthwhile project.
Thanks for the info about the augmented hessian Holger. I did notice problems with the localized orbitals when using a lot of diffuse functions, so that might be another good extension.
I added DLPNO-MP2 to the CBS driver, as well as a CBS DLPNO-MP2 test |
There was a problem hiding this comment.
Alright, a few last issues from me. I'll need to look at the sparse function one more time, once you've responded to my comments.
There are still questions for Lori about the appropriateness of the QC variables, especially when some of them are marked zero until they're added properly.
psi4/src/psi4/dlpno/mp2.cc
Outdated
| point_funcs[thread] = std::make_shared<BasisFunctions>(basisset_, grid.max_points(), nbf); | ||
| DOI_ij_temps[thread] = std::make_shared<Matrix>("(i,j) Differential Overlap Integrals", naocc, naocc); | ||
| DOI_iu_temps[thread] = std::make_shared<Matrix>("(i,u) Differential Overlap Integrals", naocc, nbf); |
There was a problem hiding this comment.
These look to me like more completely unnecessary pointers. Are these needed for parallelization, for some reason?
There was a problem hiding this comment.
We need nthread of these items so that each thread has one of its own. I changed the DOI_ij_temps and DOI_iu_temps to use regular Matrix objects. make a BasisFunctions object without a shared pointer causes a strange compilation error, so I left it as is
| * | ||
| * For all a, every value in A_to_y[a] is included in yA if at least one is present in y | ||
| */ | ||
| std::vector<int> contract_lists(const std::vector<int> &y, const std::vector<std::vector<int>> &A_to_y) { |
There was a problem hiding this comment.
Please update the documentation to include this detail, if you haven't already.
Co-authored-by: Jonathon Misiewicz <[email protected]>
doc/sphinxman/source/dlpnomp2.rst
Outdated
|
|
||
| The steep polynomial scaling (in both time and memory) of post-HF dynamic | ||
| correlation methods prohibits calculations on large systems, even for efficient | ||
| codes like |Psifour|'s :ref:`DF-MP2 <sec:dfmp2>`. This poor scaling is in part |
There was a problem hiding this comment.
| codes like |Psifour|'s :ref:`DF-MP2 <sec:dfmp2>`. This poor scaling is in part | |
| codes like |PSIfours| :ref:`DF-MP2 <sec:dfmp2>`. This poor scaling is in part |
any other lowercase |Psifour| will need to be |PSIfour|. I can't even recall how we ended up with an abbreviation harder to type than the target.
| * DLPNO-MP2 is not a drop-in replacement for DF-MP2. Instead, it should be used for | ||
| large calculations that cannot reasonably be performed with DF-MP2. The crossover | ||
| point between DF-MP2 and DLPNO-MP2 depends on details of both the calculation and | ||
| the hardware, but can be as low as 2,000 basis functions. |
There was a problem hiding this comment.
Are there restrictions on scf_type and mp2_type values, and does the program stop nicely if a forbidden one supplied? We can probably set stdsuite up to probe this in another PR if preferred.
There was a problem hiding this comment.
There are no restrictions on scf_type. DLPNO-MP2 is run density-fitted regardless of the value of mp2_type. Would we want to throw an error (or print a warning) if a user tries to run DLPNO-MP2 with mp2_type conv?
There was a problem hiding this comment.
yes, please, throw an error if mp2_type has_changed and equals CONV. better to thwart the user than surprise him. :-)
psi4/src/psi4/dlpno/mp2.cc
Outdated
|
|
||
| set_scalar_variable("MP2 SINGLES ENERGY", 0.0); | ||
| set_scalar_variable("MP2 DOUBLES ENERGY", e_mp2_corr); | ||
| set_scalar_variable("MP2 SAME-SPIN CORRELATION ENERGY", 0.0); // TODO |
There was a problem hiding this comment.
follow this pattern, I think: https://github.com/psi4/psi4/blob/master/psi4/src/psi4/dfmp2/mp2.cc#L736-L750 . If dlpno doesn't decompose readily into SS & OS, then don't set those variables.
the block you have here will get them set on the wfn. then in proc.py add this block https://github.com/psi4/psi4/blob/master/psi4/driver/procrouting/proc.py#L3932-L3934 so that global P::e has the same info.
There was a problem hiding this comment.
I can't remember the details about how the various opposite-spin and same-spin quantities should be reported, so I'll leave it to the other reviewers to comment on that. This is a really important feature and you've done an outstanding job implementing, documenting, and testing it. Kudos!
|
I've spent enough time reviewing We have three approving reviews at this point, but given the scope of the PR, we should let some time pass before merging this in. If any others want to give this PR a pass before it gets merged in, please say so now. Zach, when you've responded to Lori and Andy's comments and the PR looks ready to merge, please say so explicitly. |
|
Thanks all for the thorough reviews. I've added the OS/SS decomposition of the DLPNO-MP2 correlation energy and verified that SCS-DLPNO-MP2 works correctly. As far as I know, this PR is ready to be merged. |
|
Lori seems to have approved, and nobody else has put a hold on this PR in the week since I put out the last call, so I don't think there are objections from core devs. Fix the merge conflict, and I can merge it in once tests pass. |
|
@zachglick This PR broke the docs. Docs auto-build says Please submit a new PR to fix this. |
|
Check the formatting of the dlpno section start in read_options.cc to check it matches others, even in ways that shouldn't matter. It'll be https://github.com/psi4/psi4/blob/master/doc/sphinxman/document_options_c.pl and document_psivars that are messing up. Let me know if it proves obstinate. |
Description
Adds the DLPNO-MP2 method to Psi4. Callable as
energy('dlpno-mp2').DLPNO-MP2 (or domain-based local pair natural orbital MP2) is the method described in the following paper:
Pinski, Peter, et al. "Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals." The Journal of chemical physics 143.3 (2015): 034108
DLPNO-MP2 recovers about 99.9% of the DF-MP2 correlation energy, and scales much better in both time and memory. The following graph compares this implementation with Psi4's DF-MP2 code on linear alkanes (18 core i9-10980XE, 200 GB, cc-pVTZ basis):
Todos
Questions
DLPNObe an option for theMP2_TYPEkeyword? Should DLPNO be an option forQC_MODULE?T_CUT_DOandT_CUT_PNO. Should the other options be user-facing at all? For now, I've marked them expert. Opinions from anyone with experience running local correlation methods are welcome here.Checklist
Status