Fcuantico patch 1 by fcuantico · Pull Request #7 · NWChemEx/StructureFinder

fcuantico · 2025-06-18T17:19:16Z

Is this pull request associated with an issue(s)?
Please list issues associated with this pull request, including closing words
as appropriate.

Description
Describe what this pull request will accomplish.

TODOs
For draft pull requests please include a list of what needs to be done and check
off items as you complete them.

fcuantico · 2025-06-18T17:21:59Z

I think this how is done but not sure hopefully it will work:
So I want to propose and update to :

# Copyright 2024 NWChemEx-Project
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from simde import TotalEnergy, EnergyNuclearGradientStdVectorD


def _compare_mol_and_point(mol, points):
    """This function is essentially a work around for the comparisons not being
    exposed to Python.
    """
    if mol.size() != points.size():
        return False

    for i in range(mol.size()):
        atom_i = mol.at(i)
        point_i = points.at(i)

        for j in range(3):
            if round(atom_i.coord(j),16) != round(point_i.coord(j),16):  
                return False

    return True


def unwrap_inputs(pt, inputs):
    """ Code factorization for unwrapping a module's inputs.
    
    Many of our friends expose interfaces which are analogous to high-level
    property types like TotalEnergy, AOEnergy, and EnergyNuclearGradient.
    Furthermore, most of our friends expose all of these calculations through
    a single interface. The result is that many of our modules start by having
    dispatch logic akin to "if we're doing an energy calculation, unwrap this
    way." Or "if we're doing a gradient, unwrap this other way." This function
    factors that logic out into a single function.
    """

    mol = None
    if pt.type() == TotalEnergy().type():
        mol, = pt.unwrap_inputs(inputs)
    elif pt.type() == EnergyNuclearGradientStdVectorD().type():
        mol, point = pt.unwrap_inputs(inputs)

        if not _compare_mol_and_point(mol.molecule, point):
            raise RuntimeError(
                'Derivative must be computed at molecular geometry')
    else:
        raise RuntimeError('Property type: ' + str(pt.type()) +
                           ' is not registered')

    return mol

where the round( ,16) has been used to because I am getting the error:

press enter to continue [4.423185485e-315 4.42424859e-315 4.41864936e-315 ][0.0 0.0 1.0 ]

1: ======================================================================
1: ERROR: test_optimize_BEfire (test_optimizers.test_backward_euler_fire.Test_TotalEnergyNuclearOptimization.test_optimize_BEfire)
1: ----------------------------------------------------------------------
1: Traceback (most recent call last):
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/tests/python/unit_tests/test_optimizers/test_backward_euler_fire.py", line 48, in test_optimize_BEfire
1:     egy, pts = mm.run_as(TotalEnergyNuclearOptimization(), "BackwardEulerFire",
1:                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/src/python/structurefinder/fire/backward_euler/backward_euler.py", line 120, in run_
1:     backwardeuler_FIRE(molecule,numb_coord, R_xyz)
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/src/python/structurefinder/fire/backward_euler/backward_euler.py", line 101, in backwardeuler_FIRE
1:     G_0 =  grad_func(mol_coord,mol)           #<-- Initial Gradient
1:            ^^^^^^^^^^^^^^^^^^^^^^^^
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/src/python/structurefinder/fire/backward_euler/backward_euler.py", line 69, in grad_func
1:     return submods["Gradient"].run_as(EnergyNuclearGradientStdVectorD(), chemist.ChemicalSystem(current_molecule), current_points)
1:            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/build/_deps/friendzone-src/src/python/friendzone/nwx2qcengine/__init__.py", line 95, in run_
1:     return _run_impl('gradient', inputs, self.results(),
1:            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/build/_deps/friendzone-src/src/python/friendzone/nwx2qcengine/__init__.py", line 36, in _run_impl
1:     mol = unwrap_inputs(pts[driver], inputs)
1:           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1:   File "/home/leothan/Documents/universities/ISU/PhD/Research/NWChemX/StructureFinder/build/_deps/friendzone-src/src/python/friendzone/utils/unwrap_inputs.py", line 55, in unwrap_inputs
1:     raise RuntimeError(
1: RuntimeError: Derivative must be computed at molecular geometry
1: 
1: ----------------------------------------------------------------------

and by using round( , 16) is fixed

jwaldrop107

I know this is still a work in progress, but there are some cleanliness things to deal with.

Testing/Temporary/CTestCostData.txt

Testing/Temporary/LastTest.log

tests/python/unit_tests/test_optimizers/test_pybernyop.py

jwaldrop107 · 2025-06-18T17:42:25Z

Regarding the suggested change to friendzone/utils/unwrap_inputs.py, you'll need to open a PR in that repository with the recommended changes.

fcuantico · 2025-06-18T18:05:42Z

I think I now did the change in the FriendZone repository

…er into fcuantico-patch-1

jwaldrop107

I know this is still a work in progress, but the following are some notes on needed changes that I wanted to highlight.

jwaldrop107 · 2025-07-02T23:47:54Z

src/python/structurefinder/fire/backward_euler/__init__.py

+import pluginplay as pp
+from simde import EnergyNuclearGradientStdVectorD, TotalEnergy, MoleculeFromString
+from berny import Berny, geomlib
+import chemist
+import numpy as np
+import tensorwrapper as tw
+
+
+class GeomoptViaPyberny(pp.ModuleBase):
+
+    def __init__(self):
+        pp.ModuleBase.__init__(self)
+        self.satisfies_property_type(TotalEnergy())
+        self.description("Performs PyBerny optimization")
+        self.add_submodule(TotalEnergy(), "Energy")
+        self.add_submodule(EnergyNuclearGradientStdVectorD(), "Gradient")
+        self.add_submodule(MoleculeFromString(), "StringConv")
+
+    def run_(self, inputs, submods):
+        pt = TotalEnergy()
+        mol, = pt.unwrap_inputs(inputs)
+        molecule = mol.molecule
+
+        # Convert Chemist Chemical System to XYZ
+        xyz = ""
+        xyz += (str(molecule.size()) + "\n\n")
+        for i in range(molecule.size()):
+            xyz += (molecule.at(i).name + " " + str(molecule.at(i).x) + " " +
+                    str(molecule.at(i).y) + " " + str(molecule.at(i).z) + "\n")
+
+        # Loads the geometry string into the Berny optimizer
+        # object.
+        optimizer = Berny(geomlib.loads(xyz, fmt='xyz'))
+
+        for geom in optimizer:
+
+            # Converts the "Berny" geometry object to Chemical System
+            geom2xyz = geom.dumps('xyz')
+            print('Berny Geom to XYZ value: \n' + geom2xyz + '\n')
+            lines = geom2xyz.split('\n')
+            print('Lines of geom2xyz: \n' + str(lines) + '\n')
+            mol_string = '\n'.join(lines[2:])
+            print('Lines to string: \n' + mol_string + '\n')
+            xyz2chem_mol = submods["StringConv"].run_as(
+                MoleculeFromString(), mol_string)
+            print('String conversion from xyz to chem sys: \n' +
+                  str(xyz2chem_mol.nuclei) + '\n')
+            geom = chemist.ChemicalSystem(xyz2chem_mol)
+            print('Chemical system of xyz2chem_mol: \n' +
+                  str(geom.molecule.nuclei) + '\n')
+            geom_nuclei = geom.molecule.nuclei.as_nuclei()
+            geom_points = geom_nuclei.charges.point_set
+
+            # Main optimizer operation
+            energy = submods["Energy"].run_as(TotalEnergy(), geom)
+            print('Interim energy: \n' + str(energy) + '\n')
+            gradients = submods["Gradient"].run_as(
+                EnergyNuclearGradientStdVectorD(), geom,
+                geom_points.as_point_set())
+            print('Interim gradient: \n' + str(gradients) + '\n')
+            optimizer.send((energy, gradients))
+
+        opt_geom = geom.molecule.nuclei
+        print(
+            'Resulting relaxed geometry (assigned to variable opt_geom): \n' +
+            str(opt_geom))
+        # Optimized energy is of type "float"
+        e = tw.Tensor(energy)
+        print(e)
+        rv = self.results()
+        return pt.wrap_results(rv, e)
+
+
+def load_pyberny_modules(mm):
+    mm.add_module("PyBerny", GeomoptViaPyberny())


Suggested change

import pluginplay as pp

from simde import EnergyNuclearGradientStdVectorD, TotalEnergy, MoleculeFromString

from berny import Berny, geomlib

import chemist

import numpy as np

import tensorwrapper as tw

class GeomoptViaPyberny(pp.ModuleBase):

def __init__(self):

pp.ModuleBase.__init__(self)

self.satisfies_property_type(TotalEnergy())

self.description("Performs PyBerny optimization")

self.add_submodule(TotalEnergy(), "Energy")

self.add_submodule(EnergyNuclearGradientStdVectorD(), "Gradient")

self.add_submodule(MoleculeFromString(), "StringConv")

def run_(self, inputs, submods):

pt = TotalEnergy()

mol, = pt.unwrap_inputs(inputs)

molecule = mol.molecule

# Convert Chemist Chemical System to XYZ

xyz = ""

xyz += (str(molecule.size()) + "\n\n")

for i in range(molecule.size()):

xyz += (molecule.at(i).name + " " + str(molecule.at(i).x) + " " +

str(molecule.at(i).y) + " " + str(molecule.at(i).z) + "\n")

# Loads the geometry string into the Berny optimizer

# object.

optimizer = Berny(geomlib.loads(xyz, fmt='xyz'))

for geom in optimizer:

# Converts the "Berny" geometry object to Chemical System

geom2xyz = geom.dumps('xyz')

print('Berny Geom to XYZ value: \n' + geom2xyz + '\n')

lines = geom2xyz.split('\n')

print('Lines of geom2xyz: \n' + str(lines) + '\n')

mol_string = '\n'.join(lines[2:])

print('Lines to string: \n' + mol_string + '\n')

xyz2chem_mol = submods["StringConv"].run_as(

MoleculeFromString(), mol_string)

print('String conversion from xyz to chem sys: \n' +

str(xyz2chem_mol.nuclei) + '\n')

geom = chemist.ChemicalSystem(xyz2chem_mol)

print('Chemical system of xyz2chem_mol: \n' +

str(geom.molecule.nuclei) + '\n')

geom_nuclei = geom.molecule.nuclei.as_nuclei()

geom_points = geom_nuclei.charges.point_set

# Main optimizer operation

energy = submods["Energy"].run_as(TotalEnergy(), geom)

print('Interim energy: \n' + str(energy) + '\n')

gradients = submods["Gradient"].run_as(

EnergyNuclearGradientStdVectorD(), geom,

geom_points.as_point_set())

print('Interim gradient: \n' + str(gradients) + '\n')

optimizer.send((energy, gradients))

opt_geom = geom.molecule.nuclei

print(

'Resulting relaxed geometry (assigned to variable opt_geom): \n' +

str(opt_geom))

# Optimized energy is of type "float"

e = tw.Tensor(energy)

print(e)

rv = self.results()

return pt.wrap_results(rv, e)

def load_pyberny_modules(mm):

mm.add_module("PyBerny", GeomoptViaPyberny())

The contents in this file are just copy/paste from elsewhere. Everything but the license header can be deleted.

jwaldrop107 · 2025-07-02T23:48:54Z

src/python/structurefinder/fire/backward_euler/be2_DC_NVT.py

@@ -0,0 +1,129 @@
+#!/usr/bin/env python3


I think this whole file can be removed.

jwaldrop107 · 2025-07-02T23:49:19Z

src/python/structurefinder/fire/backward_euler/be2_DC_VarT.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python3


Similar to the other one, I think this whole file can be removed.

jwaldrop107 · 2025-07-02T23:51:57Z