py-side Hessian analysis #834
Conversation
|
Just found a checklist of stuff that needs to be working by end of this project, so may as well transcribe it somewhere it'll be remembered:
|
|
With the newest commits, this still has one embarrassing bug, but all test cases pass, and it's ready for review. This is as far as I intend to dismantle findif. |
|
Ok, I've added an orienting Description to the main PR box. Requesting review by @psi-rking, particularly the |
|
I haven't been using 'asarray'. Will start. Does the rotational_symmetry_number() code account for the change in mass for isotopic substitution somewhere? FYI, the 'nmwhess' read to me as "new mass-weighted hessian". I realized subsequently it meant 'non'. I've skimmed the changes overall, but haven't scrutinized vib.py itself. |
|
Yup, I haven't traced it fully recently, but I expect here is where mass changes get caught. In practical terms, the freq-isotope[12] test cases should cover the rotational_symmetry_number with some care. You're right that mass-symmetry-breaking in Molecule wasn't working at one time, though. I should probably spice up mints5 to check this thoroughly. Or were you concerned particularly with non-Abelian PG involving Hmm, I keep reading it as "nwchem hessian". Maybe time for a name change to just "hessian" with the weightedness in the docstring. |
|
Some checks in libmints/molecule.cc do use 'is_equivalent_to', e.g. 'has_inversion'. A possible problem I see is in molecule.cc (assuming this code hasn't changed): void Molecule::set_full_point_group(double zero_tol) contains but the 'equal_but_for_row_order' is a matrix operation If there is a problem, it would appear for higher n. |
|
Ok, @psi-rking, I've expanded the Any other problems you notice? |
|
Ok, @psi4/maintainers and @psi-rking, this mammoth is ready to review. I've updated the RN at the top, and I recommend reading them over before meeting the code. Nitpick away (esp. at numpy stuff that I was still learning at the start), and feel free to question the broader (esp. vibanal) logic flow. |
psi4/driver/driver.py
Outdated
| # At stationary point? | ||
| G0 = gradient(lowername, molecule=molecule, **kwargs) | ||
| nonzero_gradient = G0.rms() > 1.e-2 # pulled out of a hat | ||
| translations_projection_sound = not core.get_option('SCF', 'EXTERN') and \ |
There was a problem hiding this comment.
I realize it's pulled out of a hat, but it seems a little loose from my (very limited) intuition. I'll do some testing and post back here soon...
psi4/driver/driver.py
Outdated
| nonzero_gradient = G0.rms() > 1.e-2 # pulled out of a hat | ||
| translations_projection_sound = not core.get_option('SCF', 'EXTERN') and \ | ||
| not core.get_option('SCF', 'PERTURB_H') | ||
| rotations_projection_sound = translations_projection_sound and not nonzero_gradient |
There was a problem hiding this comment.
I really like this rotations_projection_sound mechanism
There was a problem hiding this comment.
EFP should also disqualify projection, I think
psi4/driver/driver.py
Outdated
| wfn.frequencies().print_out() | ||
| core.thermo(wfn, wfn.frequencies()) | ||
| # Project final frequencies? | ||
| translations_projection_sound = not core.get_option('SCF', 'EXTERN') and \ |
There was a problem hiding this comment.
This seems like a repeat of line 1313. In the spirit of the DRY principle, it may be a good idea to centralize this logic; being able to add new cases to only one location would be great
There was a problem hiding this comment.
Yeah, its awkward b/c (1) the former is in the inner fn hessian and there's no easy place to return the value to the outer frequency fn and (2) they're actually controlling two different projections – in hessian, whether to include rot in SALCs if findif, and in freq, whether to project rot in the vib analysis, whether analytic or findif.
As soon as c-side findif is dead, I agree a unified location (and hat value) is desirable.
There was a problem hiding this comment.
Regarding your second point, one solution might be to have two helper functions like energy_is_rotationally_invariant and energy_is_translationally_invariant, which could be called as part of the determination of the different *_projection_sound determinations, with appropriate additions to the resulting bool for the findif / hessian projection. Currently both have the same dependence on external potentials, but something like simple periodic boundary conditions would make the energy translationally, but not rotationally, invariant. My fear is that we could introduce something like PBC and forget to add it to the list of conditions in one place or another, potentially creating a mess. I'm sure you've already thought this through, but just wanted to voice my concerns so we keep this stuff in mind through the migration to py-side analysis.
psi4/driver/driver.py
Outdated
| # Project final frequencies? | ||
| translations_projection_sound = not core.get_option('SCF', 'EXTERN') and \ | ||
| not core.get_option('SCF', 'PERTURB_H') | ||
| rotations_projection_sound = translations_projection_sound and wfn.gradient().rms() < 1.e-2 # aforementioned hat |
There was a problem hiding this comment.
Similarly, the 1.0e-2 would be better if it were a variable, so that both locations are updated consistently in the event of us changing tolerance.
psi4/driver/driver.py
Outdated
| @@ -1709,6 +1765,73 @@ def frequency(name, **kwargs): | |||
| return core.get_variable('CURRENT ENERGY') | |||
|
|
|||
|
|
|||
| def vibanal_wfn(wfn, hess=None, irrep=None, molecule=None, project_trans=True, project_rot=True): | |||
| # TODO should go back to private | |||
| import numpy as np | |||
There was a problem hiding this comment.
Wouldn't be better to import this for the whole file? After all your hard work, I don't think we need math anymore, also.
There was a problem hiding this comment.
About to push a cleanup of driver and partial qcdb imports. Numpy now always top-level.
psi4/driver/driver.py
Outdated
| # TODO should go back to private | ||
| import numpy as np | ||
|
|
||
| if hess is None: |
There was a problem hiding this comment.
This is a great idea! Being able to analyze anything passed in is really cool.
There was a problem hiding this comment.
The vibanal_wfn was a minimum distance hookup btwn the driver and the python. There's another vibanal_str in tests/python/vib_analysis that gets at the analysis from another angle. I expect these functions to evolve a lot as we see what's really needed. The isotopes test case shows Hessian reuse (finally!).
|
|
||
| if dashlvleff.lower() == 'd2p4': | ||
| returnstring = '%12.6f %12.6f %12.6f %12.6f %12.6f %6d\n' % \ | ||
| (dashcoeff['s6'], | ||
| returnstring = dashformatter.format( |
There was a problem hiding this comment.
Ok, so I agree with your assertion that the .format is cleaner than the old % syntax now.
There was a problem hiding this comment.
I think im grumpy and old:
>>> fmt = "%12.6f %12.6f"
>>> fmt % (10, 5)
' 10.000000 5.000000'
There was a problem hiding this comment.
Ah, but mixed args & kwargs spec is so much cleaner with format. I still cling to printf in cpp, but we'll just have to disagree over % vs format in py. I love it. :-)
| text += fmt.format('[u A^2]', a='I_A', b='I_B', c='I_C', *rc_moi['u A^2']) | ||
| text += fmt.format('[GHz]', a='A', b='B', c='C', *rc_moi['GHz']) | ||
| text += fmt.format('[MHz]', a='A', b='B', c='C', *rc_moi['MHz']) | ||
| text += fmt.format('[cm^-1]', a='A', b='B', c='C', *rc_moi['cm^-1']) | ||
| print(text) |
There was a problem hiding this comment.
Is this missing a redirect to the outfile, or deliberately going to stdout?
There was a problem hiding this comment.
Pretty much the whole of qcdb outside those pieces that would clutter up a psi4 calc has a serious outfile/stdout problem for which I don't have a clean solution. Yes, it's just heading to stdout.
There was a problem hiding this comment.
We can always overwrite stdout to psi4 output (I think folks would be unhappy). Or provide a output handle to the functions.
There was a problem hiding this comment.
I'd really like a logging type approach where an output resource could be global. But an output handle is probably more realistic. I've got a lot of goodie, text = function_of_worth() going now, which complicates function returns. I may be searching the wrong things, but I haven't seen a good sol'n out there for py modules.
psi4/driver/qcdb/util/mpl.py
Outdated
|
|
||
| """ | ||
| import numpy as np | ||
| from matplotlib import pyplot |
There was a problem hiding this comment.
I'm guessing this isn't called by everyday users, but it may still be worth a try / except ImportError for the not entirely ubiquitous matplotlib.
There was a problem hiding this comment.
Hopefully they won't have do_plot=True w/o considering if they have anything to plot with, but yes, good point, there's several feature that require add'l packages that I've try/excepted so good form to do mpl, too.
There was a problem hiding this comment.
I've tripped myself up a time or two by assuming that mpl is installed on a given machine, before discovering it isn't. I'm not aware of any distros that ship it by default so I think a message describing how to pull it in through conda may aid novice users who thought they'd just get a plot automagically. Other stuff that's more standard is fine to leave unguarded IMHO
| return passed | ||
|
|
||
|
|
||
| def hessian_symmetrize(hess, mol): |
There was a problem hiding this comment.
Yay! Thanks for figuring out how to do this properly, so we can have a ceremonial burning of my horribly failed attempt.
psi4/driver/qcdb/vib.py
Outdated
| text.append(' pre-proj all modes:' + str(_format_omega(pre_frequency_cm_1, 4))) | ||
|
|
||
| # project & solve | ||
| mwhess_proj = np.einsum('ij,jk,kl->il', P.T, mwhess, P) |
There was a problem hiding this comment.
This is very elegant machinery, and I like the details being provided to the user the way they are
There was a problem hiding this comment.
Hmm, careful here without np.einsum optimize machinery this could scale like N^4. Might be better write np.dot(P.T, mwhess).dot℗. This counteracts Andy, but that could actually become a bottleneck!
|
|
||
| imagfreqidx = np.where(vibonly['omega'].data.imag > vibonly['omega'].data.real)[0] | ||
| if len(imagfreqidx): | ||
| print("Warning: thermodynamics relations excluded imaginary frequencies: {}".format(omega_str[imagfreqidx])) |
There was a problem hiding this comment.
If there's exactly 1 imaginary frequency, perhaps we could think about automagically reporting the Wigner tunneling correction.
There was a problem hiding this comment.
@CDSherrill already has plans of putting someone on the Grimme low-freq rotor correction. Wigner's another good one.
There was a problem hiding this comment.
I am assuming this (Wigner tunnelling or hindered rotor treatments) more or less fell through the cracks? Or do you have anyone working on it?
There was a problem hiding this comment.
It did make it on to the wish list, so not wholly lost. No one working on it. @rainli323 wanted a project and it might be suitably compact unless you had your eye on it.
| except ImportError: | ||
| msg = "\n\tPlot not generated; matplotlib is not installed on this machine.\n\n" | ||
| print(msg) | ||
| core.print_out(msg) |
| rotations_projection_sound = (translations_projection_sound and | ||
| stationary_point) | ||
|
|
||
| return translations_projection_sound, rotations_projection_sound |
There was a problem hiding this comment.
Yea, I like this too. Single place to find this kind of information.
There was a problem hiding this comment.
Nice! This is exactly what I was thinking when nagging about a centralized routine before.
| @@ -143,6 +143,8 @@ install: | |||
| - conda info -a | |||
| - conda create -q -n p4env python=$PYTHON_VER numpy cmake ci-psi4-lt psi4-lt-mp gau2grid dftd3 gcp snsmp2 scipy -c psi4/label/dev -c psi4 | |||
| - source activate p4env | |||
| - conda install hungarian networkx -c psi4 | |||
| - conda install deepdiff -c conda-forge | |||
There was a problem hiding this comment.
Are we adding this as Psi4 requirements?
There was a problem hiding this comment.
No. hungarian and networkx are very specialized, and they're only needed if you want to align two molecules and can't guarantee that the atoms map (first atom of molA plays same structural role as first atom of molB, etc.) So those are definitely in-function imports so that run-of-the-mill calcs aren't burdened.
deepdiff you actually will hit by default at present b/c there's a temporary? internal check that the Mol.to_dict are sound. It's easy enough to get, as I added it to conda-forge, but what do you think? Require or disable internal checks and some test cases if not present?
In fact, right now you have to compile hungarian yourself for most platforms. Author would like to shove it into numpy as his is in C and reputedly far faster than the np implementation. But no one's actively working on that. Btw, @dgasmith, do you know how a module that compiles against numpy/arrayobject.h chooses its numpy? My py27 conda pkg has been supplied with a np112 numpy but ends up with np112 in its conda build string but np110 in its setup.py build string. And complains at runtime about being compiled against a diff version. There's NPY_1_7_API_VERSION warnings at compile-time if that's a clue.
There was a problem hiding this comment.
Anything we can pip install I wouldn't worry too much about; however, if we start requiring libraries from conda-forge we enforce conda on everyone. That is a step I am not too comfortable with at the moment. As a note, you can get networkx/deepdiff off pip.
If he wants to get a C hungarian into NumPy, I can help him with that.
Yea, no idea on setup/distutils compilation, I try to stay far away.
| zero_ghost_fragments=zero_ghost_fragments, | ||
| nonphysical=nonphysical, | ||
| mtol=mtol, | ||
| verbose=verbose) |
There was a problem hiding this comment.
Cool, can these take in list as well?
There was a problem hiding this comment.
Yes, all inputs are array-like. You can see it in action in isapt1.
The output molrec dictionary is of fixed type, which is ndarray for atom stuff (geom, mass) and list for fragment stuff (fragment_separators, fragment_charges). I didn't quite understand your slack comment about which is preferable for the to_json operation (not in the scope of this PR).
| # replaced by `cls.fn = qcdb.Molecule.fn` and in qcdb/molecule.py | ||
| # itself, the raw, staticmethod fns can use their official names again | ||
| # and not need the append line at bottom of file. what I do for you, | ||
| # py2 ... |
There was a problem hiding this comment.
It might be worth considering dropping Py2 before too much longer. EOL for most major projects is end of this year, not the 2020 death clock.
There was a problem hiding this comment.
Guess we can twitter/GH/forum survey to see if there's major objection to dropping py27. The static/raw workaround wasn't so ugly as I had feared (I thought I'd have to have multiple fns with the same long args list). So while I've no affection for 27, I haven't yet reached the "Shoo!" stage.
| >>> psi4.set_memory("1e5 gb") # string w/ exponent | ||
| >>> psi4.set_memory("5e5") # string w/o units | ||
| >>> psi4.set_memory(2000) # mem too small | ||
| >>> psi4.set_memory(-5e5) # negative (and too small) |
There was a problem hiding this comment.
One of these days I'll have to actually build the docs in real time and see if my docstring cleanups are actually taking affect.
psi4/driver/qcdb/align.py
Outdated
| for jat in range(nat): | ||
| alhess[iat, jat] = (self.rotation.T).dot(blocked_hess[iat, jat].dot(self.rotation)) | ||
|
|
||
| alhess = alhess[:, self.atommap] |
There was a problem hiding this comment.
Should be able to do alhess[self.atommap[:, None], self.atommap].
There was a problem hiding this comment.
Why yes, not that I knew that when I wrote it. Even cleaner, I think I can do alhess[np.ix_(self.atommap, self.atommap)]
| basisdict = qcdb.BasisSet.pyconstruct( | ||
| mol.create_psi4_string_from_molecule(), key, resolved_target, fitrole, other, return_atomlist=return_atomlist) | ||
| basisdict = qcdb.BasisSet.pyconstruct(mol.to_dict(), | ||
| key, resolved_target, fitrole, other, return_atomlist=return_atomlist) |
There was a problem hiding this comment.
Hmm, I think we can overload the __dict__ attribute so that dict(mol) will work as well.
There was a problem hiding this comment.
Would that work on psi4.core.Molecule, too, though?
There was a problem hiding this comment.
Yea, since everything is a dict all of the special Python functions like dict, +=, etc just look for a object that has the right attributes. Kinda handy.
There was a problem hiding this comment.
Still makes me nervous b/c for both Mol classes, the to_dict could be not a full picture to recreate (say, Zmat or xyz with vars). Result of to_dict is a fine, upstanding representation of the Mol and can be reversed with from_dict, but it's not necessarily the same repr of the Mol instance.
…ag-counting issues.
… go through chgmult parser.
…base updates/fixes.
|
Ok, once again rebased and comments addressed and RTG. |
|
Thank you so much for thinking of me! I gave it some thought, and don't
quite understand what is needed. Also, I think it's a bit unlikely for me
to get to it within the next two months, is that okay?
…On Wed, Apr 10, 2019 at 2:08 PM Lori A. Burns ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In psi4/driver/qcdb/vib.py
<#834 (comment)>:
> + phi_A, phi_B, phi_C = rot_const * 100 * psi_c * psi_h / psi_kb
+ q_rot = math.sqrt(math.pi) * T**1.5 / (sigma * math.sqrt(phi_A * phi_B * phi_C))
+ sm[('S', 'rot')] = 3 / 2 + math.log(q_rot)
+ sm[('Cv', 'rot')] = 3 / 2
+ sm[('Cp', 'rot')] = 3 / 2
+ sm[('E', 'rot')] = 3 / 2 * T
+ sm[('H', 'rot')] = sm[('E', 'rot')]
+
+ # vibrational
+ vibonly = filter_nonvib(vibinfo)
+ ZPE_cm_1 = 1 / 2 * np.sum(vibonly['omega'].data.real)
+ omega_str = _format_omega(vibonly['omega'].data, decimals=4)
+
+ imagfreqidx = np.where(vibonly['omega'].data.imag > vibonly['omega'].data.real)[0]
+ if len(imagfreqidx):
+ print("Warning: thermodynamics relations excluded imaginary frequencies: {}".format(omega_str[imagfreqidx]))
It did make it on to the wish list, so not wholly lost. No one working on
it. @rainli323 <https://github.com/rainli323> wanted a project and it
might be suitably compact unless you had your eye on it.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#834 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AdbCvRCQrgbenf4C_5gu4S8gOa0kR6NTks5vfiixgaJpZM4QMWN4>
.
|
Description
Includes Joe's #772, since the original point of this was to do properly the hacks needed to get the normco into the Wfn. This includes py-side vibrational analysis, thermo, and small-system reordering. Also fixes Hessian symmetrization (py-side; haven't moved it back c-side).
This isn't a final product. I only re-hooked-back-up one conventional freq test case. There's plenty of organizational, naming, storage (want both freq and normco in Wfn? and if so, only vibs, not RT? and how ID imag freq if we're float and don't want-?) choices yet. What this does have is a substantial py test case with 6 Hessian readouts from Cfour and 6 matching p4 findif-by-grad vibrational analyses to help debug some analytic Hessian errors at wwdc.And watch out – any change tofindif.hrecompiles most of psi.In the end, this added a py-side vib analysis, identified two Hessian bugs, upgraded the two Molecule classes, added Molecule serialization everywhere except initial string input parsing, added alignment and BFS tools, and lots more – see below.
Todos
Notable points that this PR has either accomplished or will accomplish.
adds normalmode displacement export to the api and optional saving of normal modes to the wavefunction.qcdb.Moleculegeometry(np_out=True)andxyz(np_out=True)are avail to skip the cast-back-to-list step. Also,to_arraysreturns geom, mass, elez, etc. as numpy arrays, rather than needing to iterate over natom.psi4.core.Moleculeunitsno longer exported as a property. has getter/setter, and the getter returns strings, notMolecule::GeometryUnitsobjects.add_atomto take label (in add'n to symbol) and mass number, so is a full fledged CoordEntry entry pointmass_numberstorage and accessor. In Psi, this is just a pass-through – nothing done with it. Stores isotope mass number if mass corresponds to valid nuclide, -1 otherwise or unknown.set_fragmentation_patternthat sets them all at once.input_units_to_augetter/setter; checks physical reasonableness.com_fixedfn to mirrororientation_fixed. In qcdb.Mol, also afix_com.create_molecule_from_stringwas doing all kinds of contortions with fr_types and efp_chg/mult to the extend that fragment* arrays were not the same length. Function to be retired shortly, so papered over.validate_and_fill_chgmultfor general solution at the boundaries of Mol class, but no good internally to psi4.Mol b/c in python.to_arrays,to_dict,BFS,B787,scramblefns shared btwn psi4 & qcdb Mol classes. When drop py27, can attach qcdb fns directly to psi4 and drop the @static raw_ intermed fns.from_arraysconstructor that takes minimal (geom & Z) to full geom spec as arrays (homogeneous natom-length, not heterogeneous per-atom) and thoroughly validates them and fills in defaults, emitting a standardized (to become MolSSI JSON) molrec dict. Plain fn produces dict.psi4.core.Molecule.from_arraysandqcdb.Molecule.from_arraysgo a step further and return a Mol.validate_and_fill_chgmultaddresses current problems (1) create_mol_from_string's problem of not allowing overall chgmult to be set, defaulting all frag to0 1, and defaulting overall to fragA, (2) physical reasonableness (chg/mult/#elec compatible) not being tested for fragments at all and for overall not until the SCF code, (3) problems like chg mult misallocation #114, and (4) overall and frag chgmult getting out of whack upon set_chg/mult, extract_subsets, not having frag chgmult editable, scf singlet/doublet defaulting. Taking given tot/frag chg/mult info and reasoning out the missing info would have been massive logic for chg/mult/#elec simultaneously, so this just codes up the rules and the defaults and throwitertools.productat it. Slows down a bit by 10 frag, but there's room for performance optimization. Ideally we'd call this at every reinterpret_coordentry, but that's not feasible for psi4 Mol. A good newdev starter project would be (for impossible reconciliations) to keep track of the least offending combo, and print out a summary of what rules it's violating.reconcile_nucleusaddresses problems (1) can't specify atoms by atomic number, (2) can't specify masses by mass number, (3) faced with A, Z, E, real/ghost, and label@C_special@12.1, all as inputs, need to make sure (a) extracting as much user info as possible, (b) not allowing contradictory info, and (c) filling in everything else from periodictable. Like chgmult, fn handles this by defining rules and candidate values, and letting itertools do the rest. Expand atom label spec so can specify isotope and atomic number e.g.,Gh(27@58.933)or2H_deut.from_dictfn takes a fully validated and defaulted molrec dict and constructs a Mol. psi4 fn is in export_mints.cc, and qcdb fn is in molecule.py.to_dictfn serializes Mol into molrec dict. This fn in molecule.py is shared by psi4 & qcdb Mol classes. For the moment (not mandatory), this dict passes again throughfrom_arraysand output is compared to make sure (1) resulting dict is pure and (2) Mol hasn't been tampered with in an unphysical way.from_arraysused internally to qcdb.Molecule for decomposing and reconstituting Mol into np arrays for fragmentation and alignment (potentially reordering)molecule {...}, it's still getting exercised a lot.from_arraysrather thanmolecule {...}for a bit (until new molecule string parser in) if splitting bonds. See note intests/isapt1. @bwb314, your SAPT consumer base may need informing about this change when it hits master.driver.frequency()after the internal hessian() callWavefunction.frequenciesmember data with vib-only real frequencies in cm^-1 so tests happysymmetrize_hessianfunctionprint_molden_vibsprint_vibs(vibinfo)function in which q/w/x normco, number of vibs per row, number of decimal places for freq, number of decimal places for normco, and normco (nat x 3) or (3*nat x 1) is all customizable by optionFD_PROJECTfor user to force inclusion/disclusion.hungarianandnetworkxadd'l packages. Handles shuffled atoms up to at least benzene trimer and more if haven't too many symmetry-equivalent atoms. Really good for perfect matches, fine but untuned for imperfect matches. Also provides functions that use that return set to properly manipulate list, vec, mol, grad, hess (hess untested). Upon request, also tests mirror-image alignment.(only methane and ammonia pass)All pass now, thanks to @andysimfrequency()does a gradient to determine if rotational projection is safe, then another flag to query external field and company to see if even translational projection is appropriate.compare_dictsandcompare_molrecscomparison functions. Former usesdeepdiffmodule and does recursive comparison (sensitive to types). Can exempt fields withforgivearg. Latter makes use of former and also allows geoms to change if warranted by fix_com/orientation (checks via aligner).TestComparisonError, notsys.exit(1).0 2in some test cases, b/c the reconciler wants physically reasonable values, but the input molecule isn't currently passing through the code that would set those defaults.Questions
Good stuff still to do
closes #940
Status