From 9f53d4cfac47459f7ac59baba5fdf5cd75536b3f Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Tue, 14 Nov 2023 15:38:43 -0500 Subject: [PATCH] Fix GFN-FF energy - redirect GFN-FF output through a hack Signed-off-by: Geoff Hutchison --- .../qtplugins/forcefield/scripts/gfnff.py | 62 ++++++++++++++++--- 1 file changed, 55 insertions(+), 7 deletions(-) diff --git a/avogadro/qtplugins/forcefield/scripts/gfnff.py b/avogadro/qtplugins/forcefield/scripts/gfnff.py index e8f5b1ae73..ddb0d44596 100644 --- a/avogadro/qtplugins/forcefield/scripts/gfnff.py +++ b/avogadro/qtplugins/forcefield/scripts/gfnff.py @@ -4,16 +4,59 @@ import argparse import json import sys +import os try: import numpy as np from xtb.libxtb import VERBOSITY_MUTED - from xtb.interface import Calculator, Param + from xtb.interface import Calculator, Param, Environment imported = True except ImportError: imported = False +# we need to redirect stdout +# this is from https://eli.thegreenplace.net/2015/redirecting-all-kinds-of-stdout-in-python/ +from contextlib import contextmanager +import ctypes +import io +import tempfile + +libc = ctypes.CDLL(None) + +@contextmanager +def stdout_redirector(stream): + # The original fd stdout points to. Usually 1 on POSIX systems. + original_stdout_fd = sys.stdout.fileno() + + def _redirect_stdout(to_fd): + """Redirect stdout to the given file descriptor.""" + # Flush the C-level buffer stdout + #libc.fflush(c_stdout) + # Flush and close sys.stdout - also closes the file descriptor (fd) + sys.stdout.close() + # Make original_stdout_fd point to the same file as to_fd + os.dup2(to_fd, original_stdout_fd) + # Create a new sys.stdout that points to the redirected fd + sys.stdout = io.TextIOWrapper(os.fdopen(original_stdout_fd, 'wb')) + + # Save a copy of the original stdout fd in saved_stdout_fd + saved_stdout_fd = os.dup(original_stdout_fd) + try: + # Create a temporary file and redirect stdout to it + tfile = tempfile.TemporaryFile(mode='w+b') + _redirect_stdout(tfile.fileno()) + # Yield to caller, then redirect stdout back to the saved fd + yield + _redirect_stdout(saved_stdout_fd) + # Copy contents of temporary file to the given stream + tfile.flush() + tfile.seek(0, io.SEEK_SET) + stream.write(tfile.read()) + finally: + tfile.close() + os.close(saved_stdout_fd) + def getMetaData(): # before we return metadata, make sure xtb is in the path @@ -57,16 +100,21 @@ def run(filename): if "totalSpinMultiplicity" in mol_cjson["properties"]: spin = mol_cjson["properties"]["totalSpinMultiplicity"] - calc = Calculator(Param.GFNFF, atoms, coordinates, + # xtb doesn't properly mute + # so we redirect stdout to a StringIO + # and then just ignore it + f = io.BytesIO() + with stdout_redirector(f): + calc = Calculator(Param.GFNFF, atoms, coordinates, charge=charge, uhf=spin) - calc.set_verbosity(VERBOSITY_MUTED) - res = calc.singlepoint() + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() # we loop forever - Avogadro will kill our process when done while True: # read new coordinates from stdin for i in range(len(atoms)): - coordinates[i] = np.fromstring(input()) + coordinates[i] = np.fromstring(input(), sep=' ') # .. convert from Angstrom to Bohr coordinates /= 0.52917721067 @@ -74,8 +122,8 @@ def run(filename): calc.update(coordinates) calc.singlepoint(res) - # first print the energy of these coordinates print("AvogadroEnergy:", res.get_energy()) # in Hartree + # times 2625.5 kJ/mol # now print the gradient # .. we don't want the "[]" in the output @@ -87,7 +135,7 @@ def run(filename): if __name__ == "__main__": - parser = argparse.ArgumentParser("GFN2 calculator") + parser = argparse.ArgumentParser("GFN-FF calculator") parser.add_argument("--display-name", action="store_true") parser.add_argument("--metadata", action="store_true") parser.add_argument("-f", "--file", nargs=1)