Minimum value from non-diagonalized Hessian is replaced

[jeherr]

Just to keep track of the progress, I'm going to paste a couple back and forth emails between me and Per-Ola on this issue. First, per @peonor

I have a very strange problem when I’m using calculate.py to get the Hessian from the attached Jaguar file. The command I’m using is:

$SCHRODINGER/run /home/kcjh508/Q2MM/q2mm-master/q2mm/calculate.py -jh jag_Acc_Fine_B3LYP-D3_LACVPs.01.in -p

I’m also using another of my older tools to read out the same data. All elements are correct, except one:

h_jag_Acc_Fine_B3LYP-D3_LACVPs_26-25 gets the value “nan”! The raw value is -1.295139E-01, it should be converted to -1.205E+03.

I can’t figure out why it happens. Using my own predecessor of Eric’s script runs without a problem. Anything you can help with?

The really strange thing? If I change the last decimal to zero in the .01.in-file, I instead get the nan error on element 35_34…

[jeherr]

This is happening in datatypes.replace_minimum(). The jaguar hessian data is collected in calculate.py and datatypes.replace_minimum() is called on this line. The 26-25 hessian value is the lowest of the hessian matrix. The value is set to be equal to "value" which is by default 1 but back in calculate.py the "invert" variable is passed in for "value" instead, which comes from the argument parser. When I run this command with your file, opts.invert is "None". I'm not sure what behavior you are wanting here so I don't know how you want to fix it, but that is what is happening.

The reason that the NaN appears for the 35_34 element if you change the last decimal of 26_25 to 0 is because that makes 35_34 your new lowest values in the hessian.

[jeherr]

Per @peonor

This explains it, but it’s actually wrong behavior. The Hessian itself should never have a value replaced. The only time we replace a value is if we first diagonalize, and then replace the most negative eigenvalue. If there is a function that looks for the lowest value in the hessian itself and replaces it, that’s an error, I can never see a reason for such a behavior.

We have two major cases. Either we read in the raw Hessian, or we read in the eigenvalues/eigenmodes. In principle, these are equivalent data. For each of these, we can chose to use it as is, modify an eigenvalue, or modify weights. In my mind, the default has always been to use the unmodified version, which is appropriate for ground state parameterization. In my original code, Hessian modification was always something done external to Q2MM. However, at some point I think Eric considered the Hessian modification to be the default behavior, and moved it into the program. It seems this has led to bugs…

Currently, most of us have been using the Eigenmode fitting, so let’s go through that first. In this case, we have two sets of choices. The first is if we should modify the most negative eigenmode. The default should be to not modify; if we do, a new positive replacement value should be specified, and I think that should be explicit in the input, to make it more visible, not hidden as default behavior. The second choice is the weights to use, and for eigenmodes, we can differentiate between diagonal and off-diagonal elements (the reference data for off-diagonals is always zero, since these are normal modes). We can also chose to set the weight for the first element to zero. Right now, this seems to be the default, which is our preferred method for transition states, but not for ground states. This should be explicit in the input again. Right now, I have to go into constants.py to change it, which means I need different codes for ground and transition states…

The older method is to use Hessian data in raw form. With this data, if we want to use it for transition states, it must go through a modification procedure. We can chose whether or not to use mass-weighting (it is the default), then we must diagonalize it (hard, my code for that frequently crashed, so I had it in external code and did it once only), identify the most negative eigenmode, change it to a large positive value, and back-transform it to Hessian format. I’m not 100% sure Eric actually understood this procedure in detail, not if he simply looked for the most negative value in the Hessian and replaced that… And for ground states, no modification should be made, so if we have the modification there, it should not be the default. Here, we also have a problem with dummy atoms, which should be removed from the Hessian since QM can’t do them right. (the dummy atom could be it’s own discussion, longer than this…)

For the Hessian, we have a few choices for the weights. One recent development was to set the weights of block diagonal elements to zero. The motivation is that the diagonal elements are simply the sums of all other elements, so they were redundant, and also led to coupling between remote atoms in the structure. My original code had another option (taken from Hagler), to use weights between Hessian elements based on the number of bonds in between them. Hessian elements always connect two atoms. If it’s the same atom, then it’s a block diagonal element, and should be removed according to the new principles. If the atoms are directly bonded, the Hessian elements define a bond force constant, so the weight should be high. If they are separated by two bonds, that’s an angle, so it could be lower. Three bonds means a torsion, and here I followed Hagler and used a much higher weight, to allow the torsions to be fit from Hessian, not from equilibrium structure. More than three bonds mean it’s a non-bonded interaction, so unless we’re letting Hessian elements influence charges, these could be set to zero. But the option isn’t there in current Q2MM.

Happy to discuss this further, but I think the current Q2MM is wrong for a few of these choices, lack some, and the defaults can be discussed. I don’t like having to go into the code and change depending on whether I want to use it for ground or transition states.

[ericchansen]

• I think some things should stay in constants.py---things that are actually constants in that they never should change (ex. atomic massses, unit conversions, etc.).

  Anything that's actually an input should go into a JSON. We should include several examples of these JSON files in an examples directory.

  We'd have to write code for reading these JSON inputs. I just never could make time for these ease of use issues.

  What should go in the JSON input is shown below.

  • Value to replace most negative eigenvalue with.
  • All weights (remembering that there are multiple weights for various eigenmodes).
  • Step sizes.

• We should remove

  q2mm/q2mm/calculate.py

  Line 1458 in 10796ba

  datatypes.replace_minimum(hess, value=invert)

  from jh. This is definitely a bug. Look at the code directly above this line.

  q2mm/q2mm/calculate.py

  Lines 1454 to 1458 in 10796ba

  	if invert:
  	evals, evecs = np.linalg.eigh(hess)
  	datatypes.replace_minimum(evals, value=invert)
  	hess = evecs.dot(np.diag(evals).dot(evecs.T))
  	datatypes.replace_minimum(hess, value=invert)

  Under if invert: is where the inversion should be taking place (and is). This section also correctly converts from the Hessian to the eigenmodes, replaces the most negative eigenmode and then reforms the Hessian.

  I originally added the bug here 582d783. I was halfway to fixing the bug here f0b7c3c. It looks like it has stayed the same ever since. My b.

• By default, jh uses the Hessian matrix without modifications. You must specify --invert to evoke the inversion behavior (well, that is once the bug I mentioned above is fixed).

• We should create separate issues for enhancements.

  • We really need to restructure inputs going to calculate.py. It's getting out of hand. This should be discussed in a separate feature issue.
  • Create a JSON input file. Issue created Add JSON input file #67.
  • Alternative weighting scheme for Hessian. Issue created Alternative weighting scheme for Hessian elements based on Hagler #68.

[peonor]

This is still a problem. If I use compare.py with "-jh" for reference data and "-mh" for calculated data, the program crashes. There seems to be two problems. In the code, the weights are taken from "c" in this case (line 165 in the code). If I change "c" to "r" in this line, it runs, but the most negative value in the QM Hessian is replaced by "nan"! I only discovered this because I forgot to include the "--invert", but it's incorrect behavior which will cause ground state parameterizations to crash.

[v3op01]

Can you share simple example files I can test my macromodel hessian weight assignment scheme?

[peonor]

The two examples I’ve been working with here are either too simple or too complex. I’ll see if I can come up with a Jaguar example here, but in the meantime, could you just ask Jess for any single structure from one of her force field optimizations, and use that to test out the MacroModel-Gaussian combination? Preferably the smallest she has from the Pd-allyl TS force field. /Per-Ola From: Kevin Koh <notifications@github.com> Sent: den 24 oktober 2019 17:40 To: Q2MM/q2mm <q2mm@noreply.github.com> Cc: Norrby, Per-Ola <Per-Ola.Norrby@astrazeneca.com>; State change <state_change@noreply.github.com> Subject: Re: [Q2MM/q2mm] Minimum value from non-diagonalized Hessian is replaced (#57) Can you share simple example files I can test my macromodel hessian weight assignment scheme? — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub<#57?email_source=notifications&email_token=ADKVENBT6IAUQ66NRCTLPA3QQG6V5A5CNFSM4HCKMGM2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECFPKQQ#issuecomment-545977666>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADKVENBOWCU5AHBCR5BDUJLQQG6V5ANCNFSM4HCKMGMQ>.

…

________________________________ 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.

[v3op01]

the problem at the moment is that c.wht has to be assigned for the different hessian matrix. For amber, I export the 12,13,14 interation info file during its calculation and imports it when Datum is building in calculate.py, which I later call in compare.py using c.wht. For Macromodel, I did not do this, which became a problem calling c.wht since it was not assigned such parameter. I will push the fix soon.

[peonor]

While you are doing that, could you also export the 11-interaction? Then you wouldn’t need to make a separate test, just setting the 11-interaction weight to zero will remove the block diagonal elements, simplifying the code. From: Kevin Koh <notifications@github.com> Sent: den 24 oktober 2019 17:50 To: Q2MM/q2mm <q2mm@noreply.github.com> Cc: Norrby, Per-Ola <Per-Ola.Norrby@astrazeneca.com>; State change <state_change@noreply.github.com> Subject: Re: [Q2MM/q2mm] Minimum value from non-diagonalized Hessian is replaced (#57) the problem at the moment is that c.wht has to be assigned for the different hessian matrix. For amber, I export the 12,13,14 interation info file during its calculation and imports it when Datum is building in calculate.py, which I later call in compare.py using c.wht. For Macromodel, I did not do this, which became a problem calling c.wht since it was not assigned such parameter. I will push the fix soon. — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub<#57?email_source=notifications&email_token=ADKVENC2MXJE6S7DR4G4DE3QQG7ZFA5CNFSM4HCKMGM2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECFQLDI#issuecomment-545981837>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADKVENFXI7BAIQN5AHSX5PTQQG7ZFANCNFSM4HCKMGMQ>.

…

________________________________ 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.