leeping / geometric Goto Github PK

I want to use Gaussian to calculate hessian matrix for the first step, and I notice geomeTRIC can read hessian.txt from disk. But what should I do to change the hessian from “Gaussian format” to “Numpy format”? Do I only need to convert the triangular matrix got from Gaussian to a symmetric matrix? Should I change the units?

pip install geometric fails, due to this error python-versioneer/python-versioneer#375

Overview

Running transition=True fails with geometric_run_json because this call (or the other similar call here) instantiate an Optimizer with dirname=None. When the calcEnergyForce calls calc_cartesian_hessian it passes self.dirname (which is None). The first lines in calc_cartesian_hessian attempt to lookup a previously calculated hessian and the os.join call fails with TypeError because it cannot join None as a valid directory path:

geomeTRIC run_json error:
Traceback (most recent call last):
  File "/opt/conda/lib/python3.9/site-packages/geometric/run_json.py", line 253, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, None, params)
  File "/opt/conda/lib/python3.9/site-packages/geometric/optimize.py", line 705, in Optimize
    return optimizer.optimizeGeometry()
  File "/opt/conda/lib/python3.9/site-packages/geometric/optimize.py", line 623, in optimizeGeometry
    self.calcEnergyForce()
  File "/opt/conda/lib/python3.9/site-packages/geometric/optimize.py", line 232, in calcEnergyForce
    self.Hx0 = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=True, verbose=self.params.verbose)
  File "/opt/conda/lib/python3.9/site-packages/geometric/normal_modes.py", line 72, in calc_cartesian_hessian
    hesstxt = os.path.join(dirname, "hessian", "hessian.txt")
  File "/opt/conda/lib/python3.9/posixpath.py", line 76, in join
    a = os.fspath(a)
TypeError: expected str, bytes or os.PathLike object, not NoneType

Will create a PR to fix this...

When dealing with geometry optimization (no matter what particular engine is in use) I often find myself wanting to change some settings on the fly. If calculation is cheap it could be done by restarting from scratch with new settings, but when each point requires 30 minutes or more it becomes a huge waste. Breaking ongoing optimization and properly restarting it is always a pain.

But what if we could intervene into running optimization? Possible features include:

• return to the previous step of selected number;
• increase or decrease trust radius;
• freeze trust radius at specific value, or let it change freely again;
• force rebuild of Hessian and/or coordinates;
• "patch" existing Hessian, for example replace negative eigenvalue with zero or positive in minimum search;
• change number of used CPU cores (threads);
• pause optimization gracefully after current gradient calculation finishes so it frees up most of used memory and could be restarted later after some other task;
• change convergence criteria or constraints;
• change logging verbosity;
• dump current Hessian to file.

I can imagine at least a couple of ways how it could be implemented:

• add interactive command line interface to geometric-optimize, allowing interactive input of commands from user controlling terminal;
• listen on unix or tcp socket and allow client to attach, internally it could accept either raw commands like in a previous approach (so user can use anything as simple as nc to enter commands), or use some kind of wire protocol and rely on custom client application to handle user interaction.

There are a few empirical corrections which are often used in DFT computations and depend solely on geometry of the system, for example DFT-Dn (with n < 4 [*]), gCP, DFT-C. Normally, they are computed by underlying QC code and included into values of energy and gradient. However, I think it might be beneficial to treat correction value (or sum of correction values) separately from electronic energy when computing steps because of the following:

• near minimum changes in values of electronic energy and empirical correction very often became comparable and go into opposite directions;
• electronic energy and its gradient may be affected by multiple sources of numerical noise, e.g. imperfect SCF convergence, low precision of integral calculations, numeric approximations like COSX, etc., while at the same time empirical corrections and their gradients are defined by analytical expressions and are always precise;
• empirical corrections are very cheap as compared to DFT and can be easily computed for target geometry before step is made.

For example I guess it shouldn't be too expensive to calculate empirical correction for target geometry on each iteration when step is optimized to fit into trust radius. Another approach I can imagine is using a separate Hessian for correction, which could be exactly recomputed whenever it's needed.

Alternatively, this feature could be helpful simply for using empirical corrections which are unavailable in particular QC code. For example, DFT-C is not available in most QC packages I'm aware of, and there are packages that don't implement gCP and/or DFT-D4.

[*] DFT-D4 depends on partial atomic charges so its case is more complicated

Hi,

I am using geomeTRIC in the following way
geometric-optimize --psi4 input.in

How do I ensure that I am using all the cores of my system? Are there any ways to reduce the duration it takes to optimize the structures?

Regards,
Hemant

For some QCFractal OptimizationDataset workflows, it would be very useful to take at most K optimization steps and still return successful completion. For example, when starting from snapshots sampled from an MM surrogate potential at 300K, the initial QM energies may be too high, but a few steps of optimization is sufficient to explore the thermodynamically relevant energy surface in the vicinity of the starting geometry.

QCFractal / QCEngine appears to call geomeTRIC in a manner that passes in a JSON document via run_json. Options such as whether to use tric or cart appear to be passed this way. If an option such as success_on_maxiter could be passed to signal that reaching maxiter should still successfully return completion of the optimization trajectory, this use case could easily be supported.

As discussed 3/12/2019 with @leeping, it is useful to implement a consistent error code, for failures of QM gradient jobs.

The current error behavior is engine-dependent, which means when QM failure happens, some engine like Psi4 will trigger a RuntimeError, some may trigger other error types.

The improved error handling should return a consistent error code when such failure happens. Something like:

Exit with code 99.
QMError: "Gradient job failed"

Another custom type of error could be something like:

Exit with code 98.
NotConvergedError: "geomeTRIC optimization failed to converge in 50 iterations"

With the improved error handling, programs that call geomeTRIC by CLI, such as torsiondrive, will be able to get information about errors happened in geomeTRIC and determine what to do next.

Dear developers,

Is there any play to the famous ORCA? A very powerful QC package.

I've recently worked with CP2K, and its optimizer is far from being stellar. AFAIU, they use fractional coordinates (XYZ projected onto lattice vectors) which is clearly suboptimal for crystals composed of organic molecules. Observable behavior is that convergence is very slo, e.g more than 200 steps were needed to optimize cell of 80 atoms.

I wonder what kind of changes would be needed to make geomeTRIC work for this case. It seems to me that TRIC coordinates should be an ideal fit for molecular crystal. However it might be needed to take into account that certain groups of atoms on the opposite sides of unit cell are in fact parts of the same molecule dissected by cell borders.

Another thing to consider is optimization of unit cell parameters which should ideally be performed in conjunction with geometry optimization. Distances between separate molecules and cell parameters are interdependent. CP2K uses stress tensor to optimize cell parameters.

When calling geomeTRIC from pyscf, I meet the following TypeError. I am not sure if it is caused by geomeTRIC itself, or pyscf calls geomeTRIC in an improper way. @sunqm

  File "/root/pyscf/pyscf/geomopt/geometric_solver.py", line 171, in kernel
    maxsteps=self.max_cycle, **self.params)
  File "/root/pyscf/pyscf/geomopt/geometric_solver.py", line 119, in kernel
    constraints=constraints, **kwargs)
  File "/root/anaconda3/lib/python3.7/site-packages/geometric-0.9.5.1+33.g0007fbf-py3.7.egg/geometric/optimize.py", line 1755, in run_optimizer
    progress = Optimize(coords, M, IC, engine, dirname, params, xyzout)
  File "/root/anaconda3/lib/python3.7/site-packages/geometric-0.9.5.1+33.g0007fbf-py3.7.egg/geometric/optimize.py", line 1312, in Optimize
    return optimizer.optimizeGeometry()
  File "/root/anaconda3/lib/python3.7/site-packages/geometric-0.9.5.1+33.g0007fbf-py3.7.egg/geometric/optimize.py", line 1271, in optimizeGeometry
    self.calcEnergyForce()
  File "/root/anaconda3/lib/python3.7/site-packages/geometric-0.9.5.1+33.g0007fbf-py3.7.egg/geometric/optimize.py", line 983, in calcEnergyForce
    self.E = spcalc['energy']
TypeError: tuple indices must be integers or slices, not str

The bug appears at 0007fbf, but everything works well if I revert to ed832c5.

My pyscf is at pyscf/pyscf@98a0bca.

I'm having trouble with a calculation - geometric dies with this message: Optimizer.optimizeGeometry() failed to converge. (This has only happened for this specific molecule; in general optimization has been very robust.)

Here's the trace from AWS CloudWatch. This is at the B97-D3/def2-TZVP level of theory. The QM engine was only called 18 times, for the first 18 updates, and then geometric just keeps taking steps until it crashes (there are some divide by zero errors below). The optimization seems done, both visually and in terms of ∆E/gradient after 18 steps.

(We're using geometric.internal.DelocalizedInternalCoordinates with connect=True, although this is a six-atom one-molecule system so I doubt that matters much?)

Any advice? Don't want to share calculation publicly, but happy to discuss more offline.

Hi,

I'm new with geomeTric. I'm trying to use geometric with a system that contains virtual site in OpenMM, but I get this error:
Traceback (most recent call last): File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1878, in main run_optimizer(**vars(args)) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1808, in run_optimizer progress = Optimize(coords, M, IC, engine, dirname, params) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1330, in Optimize optimizer = Optimizer(coords, molecule, IC, engine, dirname, params) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 908, in __init__ self.H0 = self.IC.guess_hessian(self.coords) File "/usr/local/lib/python3.7/site-packages/geometric/internal.py", line 3244, in guess_hessian Hprim = self.Prims.guess_hessian(coords) File "/usr/local/lib/python3.7/site-packages/geometric/internal.py", line 2521, in guess_hessian elem1 = min(Elements.index(self.elem[ic.a]), Elements.index(self.elem[ic.b])) ValueError: 'M' is not in list Traceback (most recent call last): File "/usr/local/bin/geometric-optimize", line 10, in <module> sys.exit(main()) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1878, in main run_optimizer(**vars(args)) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1808, in run_optimizer progress = Optimize(coords, M, IC, engine, dirname, params) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 1330, in Optimize optimizer = Optimizer(coords, molecule, IC, engine, dirname, params) File "/usr/local/lib/python3.7/site-packages/geometric/optimize.py", line 908, in __init__ self.H0 = self.IC.guess_hessian(self.coords) File "/usr/local/lib/python3.7/site-packages/geometric/internal.py", line 3244, in guess_hessian Hprim = self.Prims.guess_hessian(coords) File "/usr/local/lib/python3.7/site-packages/geometric/internal.py", line 2521, in guess_hessian elem1 = min(Elements.index(self.elem[ic.a]), Elements.index(self.elem[ic.b])) ValueError: 'M' is not in list.

I save the system like this from openmm:
with open('system.xml', 'w') as outfile: outfile.write(XmlSerializer.serialize(system))
and I use the same pdb file. Also I add a constraints file, because my molecules are rigid.

Do you have any advice?

Thank you in advance

Anthoni

Currently if OpenMM cannot be found a logger.info is thrown which may confuse some users as it does not harm canonical operations. Can we either change this to logger.warning, eliminate the warning, or similar instead which would be more consistent with other engines (e.g., cannot find Molpro). Line is here.

Hi,

Thanks for the package. Great stuff.

Wanted to draw your attention to a pip install issue with python 3.12 and pip 23.3.1. Got the following upon running pip install geometric:

pip install geometric
Collecting geometric
  Downloading geometric-1.0.1.tar.gz (429 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 429.1/429.1 kB 1.8 MB/s eta 0:00:00
  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  error: subprocess-exited-with-error

  × Getting requirements to build wheel did not run successfully.
  │ exit code: 1
  ╰─> [31 lines of output]
      /tmp/pip-install-3k1jekti/geometric_9eb4beb194884d6fb49a080a066f467c/versioneer.py:421: SyntaxWarning: invalid escape sequence '\s'
        LONG_VERSION_PY['git'] = '''
      Traceback (most recent call last):
        File "/home/lucaami/dev/allVenvs/comp_chem/lib/python3.12/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 353, in <module>
          main()
        File "/home/lucaami/dev/allVenvs/comp_chem/lib/python3.12/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 335, in main
          json_out['return_val'] = hook(**hook_input['kwargs'])
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/home/lucaami/dev/allVenvs/comp_chem/lib/python3.12/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 118, in get_requires_for_build_wheel
          return hook(config_settings)
                 ^^^^^^^^^^^^^^^^^^^^^
        File "/tmp/pip-build-env-6h83n5py/overlay/lib/python3.12/site-packages/setuptools/build_meta.py", line 355, in get_requires_for_build_wheel
          return self._get_build_requires(config_settings, requirements=['wheel'])
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/tmp/pip-build-env-6h83n5py/overlay/lib/python3.12/site-packages/setuptools/build_meta.py", line 325, in _get_build_requires
          self.run_setup()
        File "/tmp/pip-build-env-6h83n5py/overlay/lib/python3.12/site-packages/setuptools/build_meta.py", line 507, in run_setup
          super(_BuildMetaLegacyBackend, self).run_setup(setup_script=setup_script)
        File "/tmp/pip-build-env-6h83n5py/overlay/lib/python3.12/site-packages/setuptools/build_meta.py", line 341, in run_setup
          exec(code, locals())
        File "<string>", line 30, in <module>
        File "/tmp/pip-install-3k1jekti/geometric_9eb4beb194884d6fb49a080a066f467c/versioneer.py", line 1480, in get_version
          return get_versions()["version"]
                 ^^^^^^^^^^^^^^
        File "/tmp/pip-install-3k1jekti/geometric_9eb4beb194884d6fb49a080a066f467c/versioneer.py", line 1412, in get_versions
          cfg = get_config_from_root(root)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/tmp/pip-install-3k1jekti/geometric_9eb4beb194884d6fb49a080a066f467c/versioneer.py", line 342, in get_config_from_root
          parser = configparser.SafeConfigParser()
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      AttributeError: module 'configparser' has no attribute 'SafeConfigParser'. Did you mean: 'RawConfigParser'?
      [end of output]

  note: This error originates from a subprocess, and is likely not a problem with pip.
error: subprocess-exited-with-error

× Getting requirements to build wheel did not run successfully.
│ exit code: 1
╰─> See above for output.

Managed to successfully install from source so no problem but simply an FYI.

Perhaps it's a versioning issue with PyPI as I noticed SafeConfigParser() is not in the latest version and #174 was only merged 3 weeks ago whereas the latest released version on PyPI (1.0.1) is from May 2023.

I tried to install the geometric package using pip install geometric through pyscf[geomopt], but got the following error message.

Collecting geometric>=0.9.7.2 (from pyscf[geomopt])
  Using cached geometric-1.0.1.tar.gz (429 kB)
  Preparing metadata (setup.py) ... error
  error: subprocess-exited-with-error
  
  × python setup.py egg_info did not run successfully.
  │ exit code: 1
  ╰─> [18 lines of output]
      /private/var/folders/s3/3_6505pn269_dtgjnx_5zxxr0000gn/T/pip-install-4l9o2yyd/geometric_72cb9561ac8d4622af2dba861b9a0dcf/versioneer.py:421: SyntaxWarning: invalid escape sequence '\s'
        LONG_VERSION_PY['git'] = '''
      Traceback (most recent call last):
        File "<string>", line 2, in <module>
        File "<pip-setuptools-caller>", line 34, in <module>
        File "/private/var/folders/s3/3_6505pn269_dtgjnx_5zxxr0000gn/T/pip-install-4l9o2yyd/geometric_72cb9561ac8d4622af2dba861b9a0dcf/setup.py", line 30, in <module>
          version=versioneer.get_version(),
                  ^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/s3/3_6505pn269_dtgjnx_5zxxr0000gn/T/pip-install-4l9o2yyd/geometric_72cb9561ac8d4622af2dba861b9a0dcf/versioneer.py", line 1480, in get_version
          return get_versions()["version"]
                 ^^^^^^^^^^^^^^
        File "/private/var/folders/s3/3_6505pn269_dtgjnx_5zxxr0000gn/T/pip-install-4l9o2yyd/geometric_72cb9561ac8d4622af2dba861b9a0dcf/versioneer.py", line 1412, in get_versions
          cfg = get_config_from_root(root)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/s3/3_6505pn269_dtgjnx_5zxxr0000gn/T/pip-install-4l9o2yyd/geometric_72cb9561ac8d4622af2dba861b9a0dcf/versioneer.py", line 342, in get_config_from_root
          parser = configparser.SafeConfigParser()
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      AttributeError: module 'configparser' has no attribute 'SafeConfigParser'. Did you mean: 'RawConfigParser'?
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
error: metadata-generation-failed

× Encountered error while generating package metadata.
╰─> See above for output.

note: This is an issue with the package mentioned above, not pip.
hint: See above for details.

I originally posted about the issue here and resolved it by downgrading to Python 3.11.

I cloned the repository and then I did git pull to get the latest code. With this latest code, I am getting the following error:
Traceback (most recent call last):
File "/home/Prasanta/anaconda3/bin/geometric-optimize", line 33, in
sys.exit(load_entry_point('geometric==1.0.1+226.g72443c8.dirty', 'console_scripts', 'geometric-optimize')())
File "/home/Prasanta/anaconda3/lib/python3.8/site-packages/geometric-1.0.1+226.g72443c8.dirty-py3.8.egg/geometric/optimize.py", line 1356, in main
File "/home/Prasanta/anaconda3/lib/python3.8/site-packages/geometric-1.0.1+226.g72443c8.dirty-py3.8.egg/geometric/optimize.py", line 1150, in run_optimizer
File "/home/Prasanta/anaconda3/lib/python3.8/logging/config.py", line 71, in fileConfig
formatters = _create_formatters(cp)
File "/home/Prasanta/anaconda3/lib/python3.8/logging/config.py", line 104, in _create_formatters
flist = cp["formatters"]["keys"]
File "/home/Prasanta/anaconda3/lib/python3.8/configparser.py", line 960, in getitem
raise KeyError(key)
KeyError: 'formatters'

how to solve this?

run_json was changed so that both logging.config.fileConfig and a new log handler are created causing geometric to dump its logging to both stdout and the string stream. I don't know why it was changed this way, but it clouds up many downstream logging files. Would it be possible to remove one of these logging handler initializations?

Simply commenting out run_json:157-159 is sufficient and is slightly concerning as the logging state is reset outside of a non main initialization function.

When optimizing with coordsys="prim", it will occasionally raise the error 'PrimitiveInternalCoordinates' object has no attribute 'conmethod' from the following line

And it seems, in prim, the conmethod was never be used

in which case the creation of PrimitiveInternalCoordinates at first time would be fine, while re-creation of it fails.

I wonder if it is the problem in optmize.py which could be fixed easily, or the so called conmethod is simply not implemented for PrimitiveInternalCoordinates yet?

Currently if there is an issue with an execution the use of os.exit is used which will exit out of all python layers, see here. This currently crashes entire QCArchive pipelines making use of geomeTRIC. The use of os.exit outside of command line tools is a very bad idea for exactly this reason.

An implementation that would give you the same functionality here:

try:
    run_optimizer(**vars(args))
except MyError:
    os.exit(51)
...

geomeTRIC will fail to converge if the underlying code is computing the incorrect gradient. The failure mode is that the optimizer takes very small, bad-quality steps.

When this happens, geomeTRIC should trigger a finite-difference test of the Cartesian gradient. This will have the benefit of aiding troubleshooting, i.e. if the Cartesian gradient is wrong, then it shows the underlying gradient code as the source of the problem, and not geomeTRIC itself.

I would like to run a minimization with al heavy atoms fixed (only H atoms optimized).
Is there an example of how to do this with geometric?
Thanks,
Alberto

I'm trying to run geometric with the JSON interface. Here is the input JSON:

 qc_schema_input = {
  7         "schema_name": "QC_JSON",
  8         "schema_version": 0,
  9         "driver": "gradient",
 10         "molecule": {
 11             "geometry": [
 12                 0.31914281845092773,
 13                 -1.093637466430664,
 14                 -1.5644147396087646,
 15                 0.09283685684204102,
 16                -0.7512494325637817,
 17                 -0.10052239894866943,
 18                 -0.09279406070709229,
 19                 0.7513599395751953,
 20                 0.1004934310913086,
 21                 -0.290915846824646,
 22                 1.0967583656311035,
 23                 1.5677818059921265,
 24                 -0.5198796391487122,
 25                 -0.7551677227020264,
 26                 -2.180879592895508,
 27                 1.2282583713531494,
 28                -0.6067847013473511,
 29                 -1.931321620941162,
 30                 0.43323153257369995,
 31                 -2.1737723350524902,
 32                 -1.7016046047210693,
 33                 0.9486593008041382,
 34                 -1.1027858257293701,
 35                 0.48745453357696533,
 36                 -0.7917245626449585,
 37                 -1.287994146347046,
 38                 0.26173627376556396,
 39                 -0.960730254650116,
 40                 1.099387764930725,
 41                 -0.47136321663856506,
 42                 0.7833671569824219,
 43                 1.2891098260879517,
 44                 -0.28042128682136536,
 45                 -0.42176103591918945,
 46                 2.176612377166748,
 47                 1.6888495683670044,
 48                 0.5750753879547119,
 49                 0.7906937599182129,
 50                 2.16329026222229,
 51                 -1.1788761615753174,
 52                 0.5997258424758911,
 53                 1.971407413482666
 54             ],  
 55             "symbols": [
 56                 "C",
 57                 "C",
 58                 "C",
 59                 "C",
 60                 "H",
 61                 "H",
 62                 "H",
 63                 "H",
 64                 "H",
 65                 "H",
 66                 "H",
 67                 "H",
 68                 "H",
 69                 "H"
 70             ]
 71         },
 72         "driver": "gradient",
 73         "model": {"method": "hf", "basis": "sto-3g"},
 74         "keywords": {},
 75         "program": "rdkit"
 76     }
 77 
 78 in_json_dict = {
 79         "schema_name": "qc_schema_optimization_input",
 80         "schema_version": 1,
 81         "keywords": {
 82             "coordsys": "tric",
 83             "maxiter": 100
 84         },
 85         "input_specification": qc_schema_input
 86     }

When I specify the psi4 program I get the following error:

Traceback (most recent call last):
  File "run_geometric_json.py", line 88, in <module>
    out_json_dict_2 = geometric.run_json.geometric_run_json(in_json_dict_2)
  File "/home/chayas/anaconda3/envs/psi4-1.2rc2/lib/python3.6/site-packages/geometric-0.8.2+42.g5f5690c-py3.6.egg/geometric/run_json.py", line 102, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, dirname, params)
  File "/home/chayas/anaconda3/envs/psi4-1.2rc2/lib/python3.6/site-packages/geometric-0.8.2+42.g5f5690c-py3.6.egg/geometric/optimize.py", line 886, in Optimize
    E, gradx = engine.calc(coords, dirname)
  File "/home/chayas/anaconda3/envs/psi4-1.2rc2/lib/python3.6/site-packages/geometric-0.8.2+42.g5f5690c-py3.6.egg/geometric/engine.py", line 703, in calc
    return self.calc_new(coords, dirname)
  File "/home/chayas/anaconda3/envs/psi4-1.2rc2/lib/python3.6/site-packages/geometric-0.8.2+42.g5f5690c-py3.6.egg/geometric/engine.py", line 698, in calc_new
    gradient = np.array(ret["return_value"])
KeyError: 'return_value'

This is what the JSON output looks like:

{'schema_name': 'QC_JSON', 
'schema_version': 0, 
'molecule': {
    'geometry': [0.31914281845092773, -1.093637466430664, -1.5644147396087646, 0.09283685684204102, -0.7512494325637817, -0.10052239894866943, -0.09279406070709229, 0.7513599395751953, 0.1004934310913086, -0.290915846824646, 1.0967583656311035, 1.5677818059921265, -0.5198796391487122, -0.7551677227020264, -2.180879592895508, 1.2282583713531494, -0.6067847013473511, -1.9313216209411623, 0.43323153257369995, -2.1737723350524902, -1.7016046047210693, 0.9486593008041382, -1.1027858257293701, 0.48745453357696533, -0.7917245626449585, -1.287994146347046, 0.26173627376556396, -0.960730254650116, 1.099387764930725, -0.47136321663856506, 0.7833671569824219, 1.2891098260879517, -0.28042128682136536, -0.42176103591918945, 2.176612377166748, 1.6888495683670044, 0.5750753879547119, 0.7906937599182129, 2.16329026222229, -1.1788761615753174, 0.5997258424758911, 1.971407413482666], 
    'symbols': ['C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H'], 
    'fix_orientation': True, 
    'fix_com': True}, 
'driver': 'gradient', 
'model': {
    'method': 'hf', 
    'basis': 'sto-3g'}, 
'keywords': {}, 
'program': 'psi4', 
'memory': 1932735283, 
'success': True, 
'raw_output': None, 
'provenance': {
    'creator': 'Psi4', 
    'version': '1.2rc2', 
    'routine': 'psi4.json.run_json', 
    'cpu': 'Intel(R) Xeon(R) CPU E5-2697 v4 @ 2.30GHz',
    'hostname': 'ls01', 
    'username': 'chayas'}, 
'psi4:qcvars': {
    'CURRENT DIPOLE X': -0.00589204557899569, 
    'CURRENT DIPOLE Y': 0.011250413621672823, 
    'CURRENT DIPOLE Z': -0.0008021024951966419,
    'CURRENT ENERGY': -144.34014427128707,  
    'CURRENT REFERENCE ENERGY': -144.34014427128707, 
    'HF TOTAL ENERGY': -144.3401442712871, 
    'NUCLEAR REPULSION ENERGY': 247.8811587307185, 
    'ONE-ELECTRON ENERGY': -664.7301173229439, 
    'PCM POLARIZATION ENERGY': 0.0, 
    'SCF DIPOLE X': -0.00589204557899569, 
    'SCF DIPOLE Y': 0.011250413621672823, 
    'SCF DIPOLE Z': -0.0008021024951966419, 
    'SCF ITERATION ENERGY': -144.34014427128707, 
    'SCF ITERATIONS': 7.0, 
    'SCF N ITERS': 7.0, 
    'SCF TOTAL ENERGY': -144.34014427128707, 
    'TWO-ELECTRON ENERGY': 272.50881432093826}, 
'return_result': [-0.3060768484680957, 0.7234438982361766, 1.9318407036014387, 0.03981452663663909, 1.342284843172246, -1.3307353583014359, -0.0008277640816439502, -1.3318412477589523, 1.3322453383233177, 0.2897062618702261, -0.6945841005920805, -1.938405625452312, 1.6460046024922417, -0.638461302007268, 1.3291451533205638, -1.82067129743742, -0.9347701663391986, 0.840745611076298, -0.2449448643818894, 2.1740910507365836, 0.3879621567637566, -1.6520729336086257, 0.7597279386868547, -1.2227520596303343, 1.7147190773785341, 1.115249416804866, -0.7891705715095088, 1.6785031200297202, -0.7507788731943492, 1.1948646854302098, -1.6945908937996765, -1.1188348411371976, 0.8223806797482305, 0.27721094637951815, -2.178341043591783, -0.3620843481047197, -1.7016269610301038, 0.5766870283586947, -1.2884448017390449, 1.7748530280204144, 0.956127398624706, -0.9075915635255932], 
'properties': {
    'calcinfo_nbasis': 30, 
    'calcinfo_nmo': 30, 
    'calcinfo_nalpha': 17, 
    'calcinfo_nbeta': 17, 
    'calcinfo_natom': 14, 
    'scf_one_electron_energy': -664.7301173229439, 
    'scf_two_electron_energy': 272.50881432093826, 
    'nuclear_repulsion_energy': 247.8811587307185, 
    'scf_dipole_moment': [-0.00589204557899569, 0.011250413621672823, -0.0008021024951966419], 
    'scf_iterations': 7, 
    'scf_total_energy': -144.34014427128707, 
    'return_energy': -144.34014427128707}, 
'wall_time': 2.1866307258605957}

When I change this line to gradient = np.array(ret["return_result"]) it works, however if I use "program": "rdkit", I get the same error but this time the KeyError is for return_results.

Any plans for supporting TURBOMOLE? For GGAs and density fitting, TURBOMOLE is probably the fastest code available...

Hi,
I discovered that "pip install geometric" and "conda install -c conda-forge geometric" produce different geometTRIC versions (0.9.7.1 and 0.9.7.2, respectively). The 0.9.7.1 version is incompatible with pySCF 1.7. No issues with 0.9.7.2.

This is the issue found in sunqm/pyscf#6

This error may not be easy to reproduce because geometric/internal.py" L3021 was not triggered in new runs

Hello, are you planning on cutting a new minor release? The current one is almost 2 years old.

Using geomeTRIC to optimise a molecule used to produce an opt.xyz file however now only a state_optim.xyz is produced. Is there a command-line argument which can be used to produce the opt.xyz again? Not sure if this is a bug or a deliberate change. Thanks.

A very small portion of constrained geometry optimizations will fail if the default settings are used. In one documented case, a large negative eigenvalue appears in the Hessian and the optimization steps become very small.

If the command line options --reset --epsilon 0.0 are used, then this pathological behavior is avoided for this case. Crank in particular should use this when calling geomeTRIC from the command line. I'm starting this issue to consider whether this should be the default option for geomeTRIC itself.

After some experimentation I've found that the following command works:

geometric-optimize --ase-class=ase.calculators.orca.ORCA --ase-kwargs="{\"orcasimpleinput\":\"B3LYP D3BJ\", \"orcablocks\":`cat blocks.inp | jq -Rsa .`}"

where blocks.inp is my part of the input file including basis sets which obviously cannot be passed inline.

Is this the right way, or something more elegant is available? Anyway it would be great to have examples of such things in documentation. I guess a lot of people just don't realize they can use their software of choice with geomeTRIC, or give up trying to connect things.

How can I prevent geomeTRIC from creating specific bonds?

For example, there is a structure in attached file: PMDA_chair_dimer2_in.xyz.txt
I'd like to break bonds N(8)—H(19) and N(27)—H(40), so that respective protons become freely floating fragments, each having 3 translational coordinates.

In other words, I'd like to get H(19) and H(40) optimized by their cartesian coordinates, and the rest of both molecules to be handled with regular TRIC coordinates.

This might seem a silly idea for this simple case, but I'm actually struggling with optimization of clusters with larger numbers of molecules of the same kind, e.g. 9 and 12. I believe that separating protons that are involved in H-bond systems would decouple their positions from rotations of "host" molecules thus making the whole coordinate system less coupled and thus improving convergence.

I've tried to achieve this with specially prepared PDB file
PMDA_chair_dimer2_in.pdb.txt and used following command

geometric-optimize --ase-class=xtb.ase.calculator.XTB --ase-kwargs='{"method":"GFN2-xTB"}' --engine ase --hessian=file:PMDA_chair_dimer2_631.hess --pdb PMDA_chair_dimer2_in.pdb PMDA_chair_dimer2_in.xyz

However, geomeTRIC seems to ignore bonding information and uses fragments 1-20 and 21-40. Note that ".txt" extensions were added by me specifically to allow GitHub to upload these files. Also, I'm using GFN2 here purely for testing reasons and intend to switch to ORCA (via ASE) for the real work.

In case that it's not currently possible I would be glad to get any pointers how could I implement this feature myself.

As far as I understand code in optimize.py, when step is rejected Hessian is not updated with gradients of the last (bad) point — there is a return statement after history rollback at https://github.com/leeping/geomeTRIC/blob/master/geometric/optimize.py#L688.

However, we computed one more point of PES, and while we take a step back it seems to me like a loss of important information if we completely throw away computed energy and gradients. Instead of only decreasing trust radius we could make a step into a better direction.

Disclaimer: I have only a vague understanding of how Hessian updates work. I guess a different algorithm may be required to make a proper update after misstep.

I noticed code for MultiDihedral restraints in geometric.internal. As I read it, this is a single-degree-of-freedom restraint for a rotatable bond, constraining the mean of the torsion angles. I can add geometric.internal.MultiDihedral restraints to a coordinate system, but I am a bit wary because this type is not included in geometric.optimize.ParseConstraints() or in any examples I could find. Is this constraint type currently supported? Thanks!

QCArchive procedures, which uses torsiondrive/geomeTRIC internally, sometimes try to drive torsions on linear geometries, perhaps from the result of incorrect filtering of molecules. The errors that result from such calculations seem to occur deeper in the codebase than might be necessary. A way to detect these problem cases earlier, with a descriptive error message, would be of great value.

These are some of the errors that appear on various molecules. I will give an example QCArchive TorsionDrive ID for each case:

TorsionDrive 18045969

Error type: unknown
Error:
geomeTRIC run_json error:
Traceback (most recent call last):
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/run_json.py", line 225, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, None, params)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1331, in Optimize
    return optimizer.optimizeGeometry()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1298, in optimizeGeometry
    self.calcEnergyForce()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1002, in calcEnergyForce
    spcalc = self.engine.calc(self.X, self.dirname)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/engine.py", line 873, in calc
    return self.calc_new(coords, dirname)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/engine.py", line 865, in calc_new
    raise QCEngineAPIEngineError("QCEngineAPI computation did not execute correctly. Message: " + ret["error"]["error_message"])
geometric.errors.QCEngineAPIEngineError: QCEngineAPI computation did not execute correctly. Message: QCEngine Unknown Error: Traceback (most recent call last):
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/qcdb/intf_dftd3/runner.py", line 131, in run_dftd3_from_arrays
    dftd3_driver(jobrec)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/qcdb/intf_dftd3/runner.py", line 183, in dftd3_driver
    dftd3_harvest(jobrec, dftd3rec)  # updates jobrec
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/qcdb/intf_dftd3/runner.py", line 301, in dftd3_harvest
    raise Dftd3Error('Unsuccessful run. Possibly -D variant not available in dftd3 version.')
psi4.driver.qcdb.exceptions.Dftd3Error: Unsuccessful run. Possibly -D variant not available in dftd3 version.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/json_wrapper.py", line 217, in run_json
    json_data = run_json_qcschema(copy.deepcopy(json_data), clean)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/json_wrapper.py", line 312, in run_json_qcschema
    val, wfn = methods_dict_[json_data["driver"]](method, **kwargs)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/driver.py", line 691, in gradient
    wfn = procedures['gradient'][lowername](lowername, molecule=molecule, **kwargs)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/procrouting/proc.py", line 2070, in run_scf_gradient
    disp_grad = ref_wfn._disp_functor.compute_gradient(ref_wfn.molecule())
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/procrouting/empirical_dispersion.py", line 236, in compute_gradient
    verbose=1)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib//python3.7/site-packages/psi4/driver/qcdb/intf_dftd3/runner.py", line 135, in run_dftd3_from_arrays
    raise RuntimeError(err) from err
RuntimeError: Unsuccessful run. Possibly -D variant not available in dftd3 version.

TorsionDrive 18886561

Error type: unknown
Error:
geomeTRIC run_json error:
Traceback (most recent call last):
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/run_json.py", line 225, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, None, params)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1331, in Optimize
    return optimizer.optimizeGeometry()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1296, in optimizeGeometry
    self.step()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1100, in step
    self.checkCoordinateSystem(recover=True, cartesian=LastForce)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 959, in checkCoordinateSystem
    raise ValueError("Cannot continue a constrained optimization; please implement constrained optimization in Cartesian coordinates")
ValueError: Cannot continue a constrained optimization; please implement constrained optimization in Cartesian coordinates

TorsionDrive 18886558

Error type: unknown
Error:
geomeTRIC run_json error:
Traceback (most recent call last):
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/run_json.py", line 225, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, None, params)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1331, in Optimize
    return optimizer.optimizeGeometry()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1296, in optimizeGeometry
    self.step()
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 1085, in step
    iopt = brent_wiki(froot.evaluate, 0.0, iopt, froot.target, cvg=0.1, verbose=params.verbose)
  File "/data/homezvol0/tgokey/.local/miniconda3/envs/qcf/lib/python3.7/site-packages/geometric/optimize.py", line 187, in brent_wiki
    raise RuntimeError('Not bracketed')
RuntimeError: Not bracketed

Versons:

geometric                 0.9.7.2                    py_0    conda-forge
psi4                      1.3.2+ecbda83    py37h06ff01c_1    psi4
qcelemental               0.14.0                     py_0    conda-forge
qcengine                  0.14.0                     py_0    conda-forge
qcfractal                 0.13.1                   py37_0    conda-forge
qcfractal-core            0.13.1                   py37_0    conda-forge

Please let me know if there is any additional information that would be helpful.

Here are two minimizations, the former performed with built-in optimizer of ORCA 5.0.4 (redundant internal coordinates, RFO) and the latter by geomeTRIC (unmodified master +
0001-Okay-patch.patch.txt). This is the same system as in #180, but high-quality Hessian (B3LYP/6-31G**) was used from start to the end to aid convergence.

defgrid3 nofinalgridx was used in both cases to improve quality of gradients as compared to default settings.

In case of geomeTRIC I've used --reset=False option, otherwise high-quality Hessian would be dropped right in the beginning because of negative eigenvalues created by Hessian update. I've also used --subfrctor 2 to avoid whole system rotations with intent to keep Hessian more relevant for the current geometry.

ORCA's inputs and outputs: PMDA_chair_dimer2_vdzp_orca.tar.gz

geomeTRIC's command.sh, inputs and outputs:
PMDA_chair_dimer2_geometric.tar.gz

It would be good if the JSON constraint input could be a list of tuples rather than a string. This would make the interface a bit more canonical JSON and lower the barrier to create this input.

Source code:

geomeTRIC/geometric/run_json.py

Line 138 in 5f2fa0b

const_rep.extend([x + 1 for x in constraint["indices"]])

Explanation:

Four constraint types are supported, namely ["xyz", "distance", "angle", "dihedral"], the current implementation is correct for the last three, such as distance 1 2 or angle 1 2 3, but is incorrect for the "xyz", because to freeze the xyz of atoms 0-2, it will generate string "xyz 1 2 3", which will trigger an error by geomeTRIC. The correct string should be "xyz 1-3" or "xyz 1,2,3"

Solution:

For the special case of constraint type "xyz", it would be useful to call the nifty.commadash() function instead of adding indices:
old:

const_rep.extend([x + 1 for x in constraint["indices"]])

new:

from geometric.nifty import commadash
...

if constraint["type"] == "xyz":
    const_rep.append(commadash(constraint["indices"]))
else:
    const_rep.extend([x + 1 for x in constraint["indices"]])

This could be a misunderstanding on my part, but are the signs of the normal modes returned by frequency_analysis important? When I run geometric.normal_modes.frequency_analysis on a Hessian matrix, I get different signs for some of the cartesian displacements depending on the machine I run it on.

For example running:

freqs, n_modes, g_tot = frequency_analysis(geom, hessian, symbols)

I'll get two different values for n_modes depending on the machine I run it on with the only (significant) difference being the sign change on the second array of values. Is this to be expected, or a bug of some kind? Thanks!

Machine 1:

[
        [
            [-1.45413605e-07, 7.31568094e-02, 3.49777695e-34],
            [-4.00280118e-01, -5.80603129e-01, 8.72597124e-50],
            [4.00282426e-01, -5.80601321e-01, -2.61779137e-49],
        ],
        [
            [-8.40792347e-07, 4.64736028e-02, -7.85031947e-33],
            [6.02403552e-01, -3.68838662e-01, -1.08372224e-49],
            [-6.02390207e-01, -3.68828204e-01, 3.25116671e-49],
        ],
        [
            [-6.99575075e-02, -7.31659905e-07, 9.73449382e-38],
            [5.55204565e-01, -4.35072838e-01, 1.17356863e-54],
            [5.55217962e-01, 4.35084451e-01, -3.52070589e-54],
        ],
    ]

Machine 2:

[
        [
            [8.91322413e-08, -7.31550007e-02, -2.14334058e-33],
            [4.00301550e-01, 5.80588425e-01, 4.08054739e-33],
            [-4.00302965e-01, 5.80587316e-01, -1.22416422e-32],
        ],
        [
            [3.81291125e-06, 4.64762049e-02, 1.98828084e-33],
            [6.02353872e-01, -3.68830371e-01, 5.75189999e-33],
            [-6.02414394e-01, -3.68877798e-01, -1.72557000e-32],
        ],
        [
            [-6.99575074e-02, 2.54028178e-06, 1.09923512e-37],
            [5.55245894e-01, -4.35098804e-01, 3.31014237e-37],
            [5.55176630e-01, 4.35058483e-01, -9.93042711e-37],
        ],
    ]

The first few lines of run_json.geometric_run_json update logging.config.fileConfig; however, if no configuration is set, i.e., no formatters are already configured, this code raises an exception. Using qcengine to execute a basic optimization results in the following:

FailedOperation(error=ComputeError(error_type='unknown_error', error_message='QCEngine Execution Error:\nTraceback (most recent call last):\n  File "/opt/conda/lib/python3.7/site-packages/qcengine/util.py", line 114, in compute_wrapper\n    yield metadata\n  File "/opt/conda/lib/python3.7/site-packages/qcengine/compute.py", line 147, in compute_procedure\n    output_data = executor.compute(input_data, config)\n  File "/opt/conda/lib/python3.7/site-packages/qcengine/procedures/geometric.py", line 44, in compute\n    output_data = geometric.run_json.geometric_run_json(input_data)\n  File "/opt/conda/lib/python3.7/site-packages/geometric/run_json.py", line 159, in geometric_run_json\n    logging.config.fileConfig(logIni,disable_existing_loggers=False)\n  File "/opt/conda/lib/python3.7/logging/config.py", line 71, in fileConfig\n    formatters = _create_formatters(cp)\n  File "/opt/conda/lib/python3.7/logging/config.py", line 104, in _create_formatters\n    flist = cp["formatters"]["keys"]\n  File "/opt/conda/lib/python3.7/configparser.py", line 958, in __getitem__\n    raise KeyError(key)\nKeyError: \'formatters\'\n'))

Essentially, the code in geomeTRIC is assuming the presence of logging formatters that may not exist.

• It should be possible to compute gradients in massively parallel way, i.e. instead of sequentially running parallel calculations for every gradient, N independent serial calculation should be run. In Perl similar thing can be done by using Parallel::ForkManager that manages process pool automatically, I'm pretty sure Python has something similar as well. Benefit is that sequential SCF calculations are more efficient, and also can take different number of iterations each — with separate processes they won't slow down each other. Another benefit is ability to use more processor cores/threads — AFAIK hybrid DFT doesn't get parallelization benefits after ~20 concurrent threads, while massively parallel approach would allow to use say 32 threads on the same machine without any synchronization overhead.
• I think a cheap alternative for the full 6N-point gradient check is possible: if we take two consequential steps from optimization and add one more point (-step) to them we can check gradient alongside direction of step. This can even be done as a part of optimization algorithm if needed. Having cheap check would allow to perform such checks casually without need to allocate resources for a large task equivalent to numerical Hessian calculation.
• Instead of printing status reports every 5 steps, some progress bar module could be used. For example, I've used tqdm, it can wrap any iterable and automatically generate nice progress bar in terminal. See e.g. annulen/vibAnalysis@334d920 for usage example.
• finite_difference_grad.py should share at least part of command line argument definitions with params.py to avoid code duplication and at the same time allow using all relevant features. For example, I had to make the following patch for using ASE engine:

diff --git a/tools/finite_difference_grad.py b/tools/finite_difference_grad.py
old mode 100644
new mode 100755
index f824f05..d51f348
--- a/tools/finite_difference_grad.py
+++ b/tools/finite_difference_grad.py
@@ -133,6 +133,16 @@ def parse_fd_args(*args):
     grp_univ.add_argument('--port', type=int, help='Work Queue port number (optional).\n')
     grp_univ.add_argument('--verbose', type=bool, default=False, help='Print extra information.\n')
 
+    grp_software = parser.add_argument_group('software', 'Options specific for certain software packages')
+    grp_software.add_argument(
+        '--ase-class',
+        type=str,
+        help='ASE calculator import path, eg. "ase.calculators.lj.LennardJones"')
+    grp_software.add_argument(
+        '--ase-kwargs',
+        type=str,
+        help='ASE calculator keyword args, as JSON dictionary, eg. {"param_filename":"path/to/file.xml"}')
+
     grp_help = parser.add_argument_group('help', 'Get help')
     grp_help.add_argument('-h', '--help', action='help', help='Show this help message and exit')
 
@@ -145,7 +155,7 @@ def parse_fd_args(*args):
 
 def main():
     args = parse_fd_args(sys.argv[1:])
-    molecule, engine = get_molecule_engine(engine=args['engine'], input=args['input'], dirname='fdgrad.tmp')
+    molecule, engine = get_molecule_engine(dirname='fdgrad.tmp', **args)
     if 'port' in args:
         createWorkQueue(args['port'], debug=False)
     coords = molecule.xyzs[0].flatten() * ang2bohr

However, other engine-specific arguments should also be handled in case users of respective engines need to check gradients. Perhaps grp_software has to be exported from params.py, and maybe some other option groups too.

On the other hand, I have experience with another module for handling command line arguments: absl.flags. It allows to define flags right in the modules where they are being used, and any script that uses those modules directly or indirectly will automatically be able to parse their flags from argv. Downside is an extra dependencies and a bit less user-friendly --help. I can explain more about it if you are interested.

I'm using drms 6e-4 dmax 1e-3 criteria for optimization (which were taken from ORCA's TightOpt), and I've got into situation where optimization got stuck in endless loop with energy slowly climbing up. Quailty was always < 0 but no steps were rejected with reasoning Not rejecting step - trust below tmin = 6.000e-05.

Unreasonably low tmin caught my attention (btw, I'm using tip of master branch). Corresponding code in params.py says that

1. For DFT tmin should actually be not smaller than 1.0e-4 because of gradient errors.
2. tmin should be smaller than Convergence_drms to avoid rejection of valid steps.

Both criteria look sensible to me. However, code that computes tmin is written as

self.tmin = kwargs.get('tmin', min(1.0e-4, self.Convergence_drms*0.1))

which explicitly contradicts criterion (1), as tmin can easily drop below 1e-4.

I think proper default value of tmin should better be chosen as

min(max(1.0e-4, self.Convergence_drms*0.1), self.Convergence_drms)

With this logic tmin may only become smaller than 1.0e-4 when Convergence_drms < 1.0e-4, but still be limited by Convergence_drms*0.1 from above for large values of Convergence_drms. For my case, it gives tmin = 1.0e-4.

I can make a PR if you think this is the right way.

Lee Ping (et al.) - this is awesome! Are you considering publishing this on PyPI? It would be great to be able to include this as a normal dependency.

As a note, NumPy is deprecating the Matrix class. Users of geometric will increasingly see code blocks with the following warnings:

/home/travis/miniconda/envs/test/lib/python3.6/site-packages/geometric/rotate.py:65: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
  ymat = np.matrix(y).T
/home/travis/miniconda/envs/test/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py:68: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
  return matrix(data, dtype=dtype, copy=False)

I'm planning to work on this after finishing the transition state feature.

Mode 1: Harmonic Energy restraint on degrees of freedom previously used for constraints

• User can provide a harmonic force constant at the end of the string, leading to a harmonic restraint being created instead of an exact constraint.
• Related feature: Scan groups. User can provide $group markers before each constraint inside of a $scan block to ensure that constraints within a group are being scanned in 1D with their values changing simultaneously. It should look like $group n_steps and the constraints within a group only need to provide initial and final values.

Mode 2: General energy modification

In geomeTRIC:

• User is required to provide a PDB file and an OpenMM system XML file --perturb perturb.xml

Why not pass a force field .xml file? Requires too much legwork to get things working. For example:

• Need to make residue templates to match PDB file
• Requires making distinction between bonded and nonbonded forces
• Not easy to assign a restraint to only one degree of freedom in the system

Build a tool in geomeTRIC that assists with creating mod.xml

• This should be an “input file” tool because interactive tools are hard to automate
• "Freeze, set, scan" are not available. Interactions that involve restraints to reference values can only use the current structure.
• The form of the restraint energy function:
  • Option 1: Restraint is harmonic and a single force constant parameter should be provided.
  • Option 2: Restraint is a custom function and all needed parameters should be provided.
    • Advantage - user can define flat bottomed potentials.
    • Disadvantage - lots of room for error.
  • Option 3: A few presets are provided (harmonic being one of them).
  • There is no option to define different parameters based on element or particle number
    (don’t want to worry about combining rules or parameter lookup tables).
• The user can create restraints for:
  • A single atom selection : Applies a harmonic restraint to selected atoms or all atoms
    • Done via CustomExternalForce.
  • A pair atom selection : Applies a harmonic restraint to all pairs from 2 selections, or all atom pairs
    • Done via CustomNonbondedForce.
  • Three atoms : Applies a harmonic restraint to selected atom triplets (angle) or all angles detected from bond graph
    • Done via CustomAngleForce.
  • Four atoms : Applies a harmonic restraint to selected atom quartets (dihedral) or all dihedrals detected from bond graph
    • Done via CustomTorsionForce.
• Custom force presets are also possible but infinite customization is only available through coding
  • Modification for hydrogen bonding prevention - this is what I actually want to do

Implementation:

• Uses OpenMM
• A new Engine should be created that contains the base Engine and an OpenMM engine to compute the energy modifier
• Optimizations should use the “base + modifier” PES but the “base” energies and forces should also be reported

Hi, I posted this question to the PySCF group and I was directed to here in hopes of solving this issue

Link: pyscf/pyscf#2031

As you see, when I run PySCF and use geomeTRIC for the optimisation with 100 maxstreps, the output from geomeTRIC does not inform me of failing to converge. I am uncertain why this is the case.

Now, I am not certain if this is related to PySCF but when I do mol_opt.converged in my code, I get the return False. Which I would say is my solution, but I have also noticed that if I do have a molecule that does converge, mol_opt.converged will still return False

Thanks for any help!

For certain molecules, we are getting this error reproducibly:
Unknown Error:
Traceback (most recent call last):
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/geometric/optimize.py", line 1893, in main
run_optimizer(**vars(args))
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/geometric/optimize.py", line 1821, in run_optimizer
progress = Optimize(coords, M, IC, engine, dirname, params)
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/geometric/optimize.py", line 1331, in Optimize
return optimizer.optimizeGeometry()
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/geometric/optimize.py", line 1299, in optimizeGeometry
self.evaluateStep()
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/geometric/optimize.py", line 1277, in evaluateStep
Eig1 = np.linalg.eigh(self.H)[0]
File "<array_function internals>", line 5, in eigh
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 1470, in eigh
w, vt = gufunc(a, signature=signature, extobj=extobj)
File "/sw-pmpv/sdk/anaconda-python/2020.07/envs/anaconda-2020/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 94, in _raise_linalgerror_eigenvalues_nonconvergence
raise LinAlgError("Eigenvalues did not converge")
numpy.linalg.LinAlgError: Eigenvalues did not converge

Unfortunately, I can't share any of the offending structures due to their confidential nature. I hope we will find a non-confidential test case soon. Anyway, do you have any hints how to fix this?

It would be a good idea for geomeTRIC to output/consume Psi4's current JSON specification found here. This will make it quite a bit easier to build input jobs and avoid parsing the Psi4 log files. A gradient example can be found [here]
(https://github.com/psi4/psi4/blob/master/tests/json/gradient/input.py).

The ultimate goal is to interface the code with the expanding MolSSI QC Schema so that it can operate with any code that follows the schema. However this schema isn't quite at the point where it can be adopted.

Hello, thanks for the lovely and clean package!

I'm curious if you've considered adding a calc_hessian call to Engine(), so that analytic Hessians might be employed where possible? I know DLFIND takes Hessian callbacks, and since most quantum chemistry programs already have default machinery for computing Hessians it might be nice to offload some of the work (unless there's some special finite difference tricks you're employing here).

corin

Dear developers,

Is there any play to the xTB? A cheap and powerful QC package.

Hi,

I am trying to use some functions in this TRIC package. Basically I would like to generate a new molecule in Cartesian coordinate given its old Cartesian coordinate and their difference in the internal coordinate. I would like to use TRIC because it gives a convenient representation (6 extra scalar numbers) of sub-molecules assuming their primitive ICs are known.

As a test, I started with a case of two water molecules (two_water.xyz). I made a translation (two_water_trans.xyz) and a translation plus rotation (two_water_trans_rot.xyz) of the second water molecule and tried to use the newCartesian function to generate the new Cartesian. The procedure is the following,

mol = Molecule(fnm='mol.xyz', ftype='xyz');
mol2 = Molecule(fnm='mol2.xyz', ftype='xyz')
tric = DelocalizedInternalCoordinates(mol, build=True, connect=False, addcart=False, conmethod=1)
tric.Prims.newCartesian(mol.xyzs[0].flatten(), tric.Prims.calcDiff(mol2.xyzs[0], mol.xyzs[0])).reshape(-1, 3)

This works for translation (two_water.xyz as mol and two_water_trans.xyz as mol2) but not for translation plus rotation (two_water.xyz as mol and two_water_trans_rot.xyz as mol2). In the latter case, the generated Cartesian coordinates break the second water molecule. I attach the three xyz files and a jupyter notebook in the zip file that could reproduce the behavior.

Any help will be appreciated!

two_water.zip