Hessian inversion

[peonor]

Where did the old technique with setting a large positive value for the negative eigenvalue go? We've mostly replaced it with the new, "natural" TSFF where we just set the weight to zero, but there are still some cases where we'd like the old technique available, and I can't find it now. Just to be clear, we do not need to create the modified Hessian, just replace the reference value for the first eigenvalue with a large positive in the reference data and use it with a non-zero weight.

The application I need it for right now is a final step where the force field has been converged to the "natural" values, taking out the force constants only for the breaking and forming bonds (and maybe some angles including them if they also ended up with too low values), and rerun ONLY those parameters with the modified Hessian. See Elaines and my 2015 article in JCC, the "F-SN2" TS, for motivation.

[peonor]

A very simple solution would be to ALWAYS set this egienvalue to the large positive value by default. We wouldn't change the "normal" operation where we try to get the natural TSFF, since the weight there is zero. Thus, to include it, we would only need to change the "eig_i" in constants.py

[ericchansen]

Rather than writing separate arguments for calculate to handle Hessian inversion (past implementation), I think we should add an option to calculate that does the inversion to all other Hessian datatypes.

Could work something like this.

python calculate.py -d somedir -gh hessian_in_archive.log -i 1.5

This would read the Gaussian Hessian as per usual, but it would change element [0, 0] to be 1.5.

Going further with this idea, the Hessian is read in many ways depending upon the filetype. However, it typically looks something like this, either stored directly on the file object or on a structure inside the file object.

hess = file_ob.hess
hess = file_ob.structures[0].hess

There is already the function datatypes.replace_minimum, which we could use. This function finds the minimum value in the array or matrix, and replaces it with whatever the argument value is. It would be absurdly simple to run this function after reading the Hessian matrix if the argument -i is given to calculate.

However, before I make that rediculously simple change, I am curious about peoples' opinions on this. Should replace_minimum actually replace the minimum value in the array/matrix, or replace element [0, 0]? That being said, this should have the same effect in all cases I'm imagining.

[ericchansen]

Also, I just realized that my changes mentioned above would prevent anyone from doing inversion on some structures, while leaving the Hessian from some other structures as they are.

Is this behavior okay with everyone, at least for now, for the sake of having this issue handled promptly (this enhancement could always be opened as a separate issue later if someone actually needs this feature)?

[ericchansen]

And lastly (since I keep thinking of things), always inverting the Hessian as you mentioned Per-Ola would cause problems for optimizing GSFFs. This would certainly problematic for one of the projects currently going on in our lab that I think you have a particular vested interest in.

[peonor]

OK, agree to the simple argument, and also, I didn't think of GSFF's, sorry... I cannot imagine a case where we only want to invert some Hessians. Only one thing: would you like to give the option of adding a second value after -i, for the weight? Every Hessian I have seen has been sorted so that [0,0] is also the smallest, so I have no strong preference, but [0,0] feels simple and robust.

[ericchansen]

Okay, attempting to make this quick change in the next 10 mins.

The additional weight option could probably be added if I thought about how argparse works. For the sake of time here, I am going to leave modifying the weight to changing the value constants.WEIGHTS['eig_i'] for now.

[ericchansen]

Now that I think about it, there are a few ways to do all of this. Could someone double check my thinking and also provide input as to what method would be ideal?

1. Read the Hessian from something, decompose it into its eigenvectors and eigenvalues, invert the eigenvalue, and reform Hessian (very easy for me to write that code).
2. Read the eigenvectors and eigenvalues from something, invert the eigenvalue, and form the Hessian. Rotations and translations have been projected out by Gaussian and Jaguar? Then I think we would have to use the QM eigenvectors and MM Hessian to determine MM eigenvalues, and then use QM eigenvectors and MM eigenvalues to form the Hessian of the appropriate size.
3. We can read the Hessian and the eigenvectors, use the eigenvectors to calculate the eigenvalues, invert the eigenvalue, and then reform the Hessian. Same idea here as before. Would use QM eigenvectors and MM Hessian to calculate MM eigenvalues, and then reform Hessian of appropriate size.

[ericchansen]

If the 1st method mentioned in my last post is okay with everyone, then I think this issue should be handled by commit f0b7c3c.

[peonor]

A couple of things:

• I think we should only work with mass-weighted Hessians. Willing to discuss this, nothing inherently better, but we have a risk of messing things up if we suddenly start mixing types.
• My prefered method for obtaining this information now is to ONLY read eigenvectors and eigenvalues from the QM log files (high precision, HPModes in Gaussian). These will always be for the mass-weighted and projected Hessian (slightly different conventions in Jaguar and Gaussian though). The Hessian itself is available, in Gaussian archive and formcheck (different orientations, dangerous) or in Jaguar 01.in files, but if we use those, we must keep track of correct orientation, do manual mass-weighting, and preferably also projection (the latter is very boring). We did this without projection in the old code, resulting in redundant information that sometimes gave unfortunate interdependencies between parameters that should not be coupled. Eliminating block-diagonals helped, but this is a better solution. I don't think that we, at any point, will need the actual QM Hessian anymore.
• From MM, we should only use the mass weighted Hessian, since this is the only thing that makes sense to combine with the eigenvectors of the mass-weighted QM Hessian. If we, at some point, decide to go for non-mass-weighted eigenvectors from QM, we need to change this, but then we need to remember doing either projection or block-diagonal elimination.
• The inversion technique does NOT require you to reform the Hessian. The only thing that is needed is to replace the first (most negative) eigenvalue with a large positive number, and change the weight from zero. Repeat, ONLY change one value and weight in the reference data set, still use the same eigenvectors and calculate the MM eigenvalues the same way.
• We might want to think of a more general mechanism where an arbitrary eigenvalue could be replaced; Goddard used a technique where he replaced all eigenvalues with the experimental numbers (from IR frequencies), and used the resulting Hessian for parameterization; this was one of my main inspirations when coming up with the inversion technique. Could be good for GS force fields if we don't trust the frequencies from DFT, but this is advanced course, no need now.

[ericchansen]

1. Agreed. Stick with one form of the energy derivatives.
2. Currently -jh and -gh read the Hessian, not the eigenvectors and eigenvalues, contrary to your preferred method. Similarly, -jeigz reads the eigenvectors and Hessian, not the eigenvalues directly. However, -geigz functions like you want, directly reading the eigenvalues. More methods should be written such that -jh, -gh and -jeigz function by only reading eigenvalues and eigenvectors. We should keep the current versions as well to compare and double check results.
3. Agreed, and that's how it is. Also, I believe the the QM eigenvectors we typically get from Jaguar and Gaussian are all mass-weighted, so no worries here.
4. For the current methods in -jh and -gh, we do need to reform the Hessian. If we instead read the eigenvectors and eigenvalues, then we would still need to form the Hessian. This datatype is, after all, direct use of the Hessian. For future reference, the weights can be set in constants.WEIGHTS, and the inverted value can be set using the -i option in calculate, which was addressed for issue #30 by commit f0b7c3ca5205301a.
5. I like the idea, but agreed we already have plenty on our plates for now. :smile:

[peonor]

When you have the eigenvectors and eigenvalues, reforming the Hessian is trivial, you can go back and forth with

H = V.W.VT W = VT.H.V

As far as I see, we don't need those methods right now, I prefer what we currently use, multiply the eigenvectors onto the MM Hessian, and compare diagonal and off-diagonal elements separately. But don't throw the old methods away, we may want them in the future... (I've changed my mind before).

[ericchansen]

Okay, so it sounds like we have to address bullet point 2 in my previous post before we can close this issue then. I will admit this isn't a high priority for me, since it's doing what we already have in different ways.

[peonor]

Agree to low priority. How about temporarily deactivating the options, with a comment on the problem, but leave the functions in the code?

[ericchansen]

I mean, is there a problem though? I don't see what we need to deactivate. Currently, we simply read the Hessian directly. I feel like that option should be left for the user.

[peonor]

Fine for me, I'm not using those options currently. As long as we don't run the risk of a random user starting to compare apples with pears.

[ericchansen]

I want to get point 2.) taken care of so we can close this issue.

There are frequencies at the end of the file. They don't have many significant figures, but they're there. I already wrote code to read these, shown here. We could simply use those eigenvalues, but these are off from the current -jeigz eigenvalues by about 4 decimal places.

Can someone else double check my units here?

Here's an example file.

X001_E1.out.txt

[peonor]

Don’t have much time now, but if you want the most number of decimals, pick the frequency, convert it to a.u. (factor in my HessMan.py program), square it, and put back the sign, should now be the eigenvalue of the mass-weighted Hessian. This works with both Gaussian and Jaguar. There is some kind of definition difference if you try to pick force constants and reduced masses from the two programs, they differ (and I don’t know the underlying reason), but the frequencies are the same for the same type of calculation.

If you pick out the force constants, you have much fewer significant digits, and you have to do some kind of program-dependent conversion.

/Per-Ola

From: Eric Hansen [mailto:notifications@github.com] Sent: den 6 juli 2016 20:40 To: ericchansen/q2mm q2mm@noreply.github.com Cc: Norrby, Per-Ola Per-Ola.Norrby@astrazeneca.com; Author author@noreply.github.com Subject: Re: [ericchansen/q2mm] Hessian inversion (#30)

I want to get point 2.) taken care of so we can close this issue.

There are frequencies at the end of the file. They don't have many significant figures, but they're there. I already wrote code to read these, shown herehttps://github.com/ericchansen/q2mm/blob/b5acb60dba4430cabadbd31b65f4489324d0a041/filetypes.py#L1016-L1018. We could simply use those eigenvalues, but these are off from the current -jeigz eigenvalues by about 4 decimal places.

Can someone else double check my units here?

Here's an example file.

X001_E1.out.txthttps://github.com/ericchansen/q2mm/files/350378/X001_E1.out.txt

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/ericchansen/q2mm/issues/30#issuecomment-230865547, or mute the threadhttps://github.com/notifications/unsubscribe/ANVSNB4Aubo7Ov50Z9uPWCjwx5iko51Hks5qS_aagaJpZM4I2BP2.

---

Confidentiality Notice: This message is private and may contain confidential and proprietary information. If you have received this message in error, please notify us and remove it from your system and note that you must not copy, distribute or take any action in reliance on it. Any unauthorized use or disclosure of the contents of this message is not permitted and may be unlawful.

[ericchansen]

Revisiting what has to be done for this issue.

1. Add method of reading Jaguar eigenvalues directly. PO suggests,

... pick the frequency, convert it to a.u. (factor in my HessMan.py program), square it, and put back the sign, should now be the eigenvalue of the mass-weighted Hessian.