VMD-L Mailing List

From: Marcelo C. R. Melo (melomcr_at_gmail.com)
Date: Mon Aug 03 2020 - 17:26:11 CDT

Dear list,
I am trying to optimize the parameters for bonds in a molecule and failing
to get reasonable values.

To isolate the issue, I tried using the simple ethanol example provided in
the FFTK tutorial, but I still get unreasonably low values for the force
constant between carbon and oxygen atoms. (relevant files for the minimal
example are in
https://drive.google.com/drive/folders/1OUUsfIjDPAY1YzuHxpYsQf6SPGuumNzb?usp=sharing
)

More specifically, after optimizing the geometry of the molecule and
calculating the Hessian (all with ORCA), I move to the "Opt. Bonded" tab
and manually enter the bonded term to be optimized (CG321 OG311), giving it
initial values close to the expected (300 kcal/mole/A**2 and 1.4 A).

The issue is that I am constantly getting values for the force constant of
around 99 kcal/mole/A**2.
The value reached in the tutorial is around 358 kcal/mole/A**2.
The relevant line in CGenFF's parameter file is:
CG321 OG311 428.00 1.4200 ! PROT methanol vib fit EMB 11/21/89

Clearly there is a large difference between the values I get and those in
the tutorial (and CGenFF). Granted, I am using ORCA and the tutorial uses
Gaussian, yet the values should be comparable.

I am using the latest VMD release and I tried both NAMD 2.13 and 2.14b.
Also tried tuning the parameters in the simulated annealing procedure, but
with no significant change (always reaching force constant values ~95).
Further optimizing the geometry also did not help.

My current suspicion is that the 2-job Hessian calculation in ORCA is not
saving the output properly, and may be throwing off the parameter fitting.
The hessian calculation (job 1) ends normally, but job 2 complains about
the geometry optimization not having converged (see file "hess.out"). I
imagine this is just a trick to write the hessian to file, but I can't see
other issues right now.

Could someone suggest a direction to debug this problem?

Best,
Marcelo