From 37311208e441971fa3c8add76940d0e2697854be Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Thu, 6 Jun 2024 10:37:45 -0700 Subject: [PATCH 1/2] BAGEL engine --- geometric/engine.py | 144 ++++++++++++++++++++++++++++++++++++++++++- geometric/errors.py | 3 + geometric/params.py | 2 +- geometric/prepare.py | 13 +++- 4 files changed, 156 insertions(+), 6 deletions(-) diff --git a/geometric/engine.py b/geometric/engine.py index db009205..008ca9d3 100644 --- a/geometric/engine.py +++ b/geometric/engine.py @@ -6,7 +6,7 @@ Authors: Lee-Ping Wang, Chenchen Song Contributors: Yudong Qiu, Daniel G. A. Smith, Sebastian Lee, Qiming Sun, -Alberto Gobbi, Josh Horton +Alberto Gobbi, Josh Horton, Nathan Yoshino Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: @@ -45,11 +45,12 @@ import numpy as np import re import os +import json from .molecule import Molecule, format_xyz_coord from .nifty import bak, au2ev, eqcgmx, fqcgmx, bohr2ang, logger, getWorkQueue, queue_up_src_dest, rootdir, copy_tree_over from .errors import EngineError, CheckCoordError, Psi4EngineError, QChemEngineError, TeraChemEngineError, \ - ConicalIntersectionEngineError, OpenMMEngineError, GromacsEngineError, MolproEngineError, QCEngineAPIEngineError, GaussianEngineError, QUICKEngineError, CFOUREngineError + ConicalIntersectionEngineError, OpenMMEngineError, GromacsEngineError, MolproEngineError, QCEngineAPIEngineError, GaussianEngineError, QUICKEngineError, CFOUREngineError, BagelEngineError from .xml_helper import read_coors_from_xml, write_coors_to_xml # Strings matching common DFT functionals @@ -1822,6 +1823,145 @@ def detect_dft(self): return True return False +class Bagel(Engine): + """ + Run a Bagel energy and gradient calculation. + """ + def __init__(self, molecule=None, exe=None, threads=None): + if molecule is None: + # create a fake molecule + molecule = Molecule() + molecule.elem = ['H'] + molecule.xyzs = [[[0,0,0]]] + super(Bagel, self).__init__(molecule) + self.threads = threads + + def load_bagel_input(self, bagel_input): + """ + Bagel .json input parser. Supports both angstroms and bohr. + """ + with open(bagel_input,'r') as bageljson: + bagel_dict = json.load(bageljson) + + elems = [] + coords = [] + readang = False + + bagel_dict['bagel'][:] = [dicts for dicts in bagel_dict['bagel'] if dicts['title'].lower() in ['molecule','force']] + + if not bagel_dict['bagel'][0]['title'].lower() in 'molecule': + raise RuntimeError("The first dictionary in Bagel json file %s must be title:'molecule'." % bagel_input) + + mole_dict = bagel_dict['bagel'][0] + force_dict = bagel_dict['bagel'][1] + + for atom_dict in mole_dict['geometry']: + elems.append(atom_dict['atom']) #extract elements + coords.append(atom_dict['xyz']) #extract coordinates + if 'angstrom' in mole_dict or 'Angstrom' in mole_dict or 'ANGSTROM' in mole_dict: + #bagel by default reads coordinates in Bohr, check to see which input is used. + try: + readang = mole_dict['angstrom'] + except: + pass + try: + readang = mole_dict['Angstrom'] + except: + pass + try: + readang = mole_dict['ANGSTROM'] + except: + pass + try: + force_dict['title'].lower() == 'force' + self.force_method = force_dict['method'][0]['title'].lower() + except: + raise RuntimeError('Bagel json file %s should have a "force" title with a specified method.' % bagel_input) + + self.M = Molecule() + self.M.elem = elems + + if readang: + self.M.xyzs = [np.array(coords, dtype=np.float64)] + elif not readang:#save coordinates in molecule object as angstroms + self.M.xyzs = [np.array(coords, dtype=np.float64)*bohr2ang] + bagel_dict['bagel'][0]['angstrom']="True" #for consistency geomeTRIC will create Bagel files that use angstrom + self.bagel_temp = bagel_dict + if self.force_method.lower in ['fci', 'zfci', 'ras', 'smith', 'relsmith']: + raise RuntimeError("Full configuration interaction and SMITH3 code are not currently supported.") + + + def calc_new(self, coords, dirname): + if not os.path.exists(dirname): + os.makedirs(dirname) + # Convert coordinates back to the xyz file + self.M.xyzs[0] = coords.reshape(-1, 3) * bohr2ang + + # Change coordinates in bagel_temp dictionary + for atom in range(self.M.xyzs[0].shape[0]): + self.bagel_temp['bagel'][0]['geometry'][atom]['xyz'] = list(self.M.xyzs[0][atom]) + # Write Json file + with open(os.path.join(dirname, 'run.json'), 'w') as outfile: + outfile.write(json.dumps(self.bagel_temp, indent=2)) + outfile.close() + try: + # Run Bagel + subprocess.check_call('export BAGEL_NUM_THREADS=%s' % self.threads, cwd=dirname, shell=True) + subprocess.check_call('export MKL_NUM_THREADS=%s' % self.threads, cwd=dirname, shell=True) + subprocess.check_call('BAGEL run.json > run.out', cwd=dirname, shell=True) + # Read energy and gradients from Bagel output + result = self.read_result(dirname) + except (OSError, IOError, RuntimeError, subprocess.CalledProcessError): + raise BagelEngineError + return result + + + def read_result(self, dirname, check_coord=None): + """ read an output file from Bagel""" + if check_coord is not None: + raise CheckCoordError("Coordinate check not implemented") + energy, gradient = None, None + bagel_out = os.path.join(dirname, 'run.out') + with open(bagel_out) as outfile: + found_grad, read_grad, found_energy, read_energy = False, False, False, False + for line in reversed(outfile.readlines()): + if '* METHOD: FORCE' in line and not found_grad: + read_grad = True + gradient = [] + elif '* Nuclear energy gradient' in line and read_grad: + read_grad, found_grad = False, True + gradient.reverse() + gradient = np.array(gradient, dtype=np.float64) + + elif read_grad and any(x in line for x in ['x','y','z']): + gradient.append(float(line.split()[1])) + + #For MP2, and HF + elif '* SCF iteration converged' in line and not (self.force_method.lower() == 'casscf' or self.force_method.lower() =='caspt2') and not found_energy: + read_energy = True + + elif read_energy and not (self.force_method == 'casscf' or self.force_method=='caspt2'): + split = line.split() + if len(split) == 4: + read_energy = False + found_energy = True + energy = float(line.split()[1]) + #For CASSCF + elif '* Second-order optimization converged' in line and self.force_method.lower() == 'casscf' and not found_energy: + read_energy = True + elif read_energy and self.force_method == 'casscf': + split = line.split() + if len(split) == 5 and split[1] == '0': + found_energy = True + read_energy = False + energy = float(split[2]) + + if energy is None: + raise RuntimeError("Bagel energy is not found in %s, please check." % bagel_out) + if gradient is None: + raise RuntimeError("Bagel gradient is not found in %s, please check." % bagel_out) + return {'energy':energy, 'gradient':gradient} + class QCEngineAPI(Engine): def __init__(self, schema, program): try: diff --git a/geometric/errors.py b/geometric/errors.py index a099f363..fbac86b4 100644 --- a/geometric/errors.py +++ b/geometric/errors.py @@ -108,6 +108,9 @@ class GaussianEngineError(EngineError): class QCEngineAPIEngineError(EngineError): pass +class BagelEngineError(EngineError): + pass + class ConicalIntersectionEngineError(EngineError): pass diff --git a/geometric/params.py b/geometric/params.py index 11397a62..6e472bc4 100644 --- a/geometric/params.py +++ b/geometric/params.py @@ -330,7 +330,7 @@ def parse_optimizer_args(*args): '"psi4" = Psi4 "openmm" = OpenMM (pass a force field or XML input file)\n' '"molpro" = Molpro "gmx" = Gromacs (pass conf.gro; requires topol.top and shot.mdp\n ' '"gaussian" = Gaussian09/16 "ase" = ASE calculator, use --ase-class/--ase-kwargs\n ' - '"quick" = QUICK\n') + '"quick" = QUICK "bagel" = Bagel\n') grp_univ.add_argument('--nt', type=int, help='Specify number of threads for running in parallel\n(for TeraChem this should be number of GPUs)') grp_jobtype = parser.add_argument_group('jobtype', 'Control the type of optimization job') diff --git a/geometric/prepare.py b/geometric/prepare.py index 2a6baa97..badfacea 100644 --- a/geometric/prepare.py +++ b/geometric/prepare.py @@ -46,7 +46,7 @@ from .ase_engine import EngineASE from .errors import EngineError, InputError from .internal import Distance, Angle, Dihedral, CartesianX, CartesianY, CartesianZ, TranslationX, TranslationY, TranslationZ, RotationA, RotationB, RotationC, CentroidDistance -from .engine import set_tcenv, load_tcin, TeraChem, ConicalIntersection, Psi4, QChem, Gromacs, Molpro, OpenMM, QCEngineAPI, Gaussian, QUICK, CFOUR +from .engine import set_tcenv, load_tcin, TeraChem, ConicalIntersection, Psi4, QChem, Gromacs, Molpro, OpenMM, QCEngineAPI, Gaussian, Bagel, QUICK, CFOUR from .molecule import Molecule, Elements from .nifty import logger, isint, uncommadash, bohr2ang, ang2bohr from .rotate import calc_fac_dfac @@ -88,7 +88,7 @@ def get_molecule_engine(**kwargs): ## MECI calculations create a custom engine that contains multiple engines. if kwargs.get('meci', None): if engine_str is not None: - if engine_str.lower() in ['psi4', 'gmx', 'molpro', 'qcengine', 'openmm', 'gaussian','quick', 'cfour']: + if engine_str.lower() in ['psi4', 'gmx', 'molpro', 'qcengine', 'openmm', 'gaussian', 'bagel', 'quick', 'cfour']: logger.warning("MECI optimizations are not tested with engines: psi4, gmx, molpro, qcengine, openmm, gaussian, quick, cfour. Be Careful!") elif customengine: logger.warning("MECI optimizations are not tested with customengine. Be Careful!") @@ -134,7 +134,7 @@ def get_molecule_engine(**kwargs): engine_str = engine_str.lower() if engine_str[:4] == 'tera': engine_str = 'tera' - implemented_engines = ('tera', 'qchem', 'psi4', 'gmx', 'molpro', 'openmm', 'qcengine', "gaussian", "ase", "quick", "cfour") + implemented_engines = ('tera', 'qchem', 'psi4', 'gmx', 'molpro', 'openmm', 'qcengine', "gaussian", "bagel", "ase", "quick", "cfour") if engine_str not in implemented_engines: raise RuntimeError("Valid values of engine are: " + ", ".join(implemented_engines)) if customengine: @@ -298,6 +298,13 @@ def get_molecule_engine(**kwargs): logger.info("The gaussian engine exe is set as %s\n" % engine.gaussian_exe) # load the template into the engine engine.load_gaussian_input(inputf) + elif engine_str == 'bagel': + logger.info("Bagel engine selected. Expecting Bagel input for gradient calculation.\n") + engine = Bagel(threads=threads) + engine.load_bagel_input(inputf) + M = engine.M + threads_enabled = True + elif engine_str == 'cfour': logger.info("CFOUR engine selected. Expecting CFOUR input for gradient calculation.\n") engine = CFOUR(inputf) From 514df1f2082abdde258d9033a1426ce5445c4e0e Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Fri, 7 Jun 2024 11:51:24 -0700 Subject: [PATCH 2/2] Add input and output files for P4 CASSCF --- .../P4_casscf_bagel/P4_casscf_ang.json | 27 ++++ .../saved/2024-06-07/P4_casscf_ang.log | 122 ++++++++++++++++++ .../saved/2024-06-07/P4_casscf_ang_optim.xyz | 30 +++++ 3 files changed, 179 insertions(+) create mode 100755 examples/1-simple-examples/P4_casscf_bagel/P4_casscf_ang.json create mode 100644 examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang.log create mode 100644 examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang_optim.xyz diff --git a/examples/1-simple-examples/P4_casscf_bagel/P4_casscf_ang.json b/examples/1-simple-examples/P4_casscf_bagel/P4_casscf_ang.json new file mode 100755 index 00000000..e15e78d3 --- /dev/null +++ b/examples/1-simple-examples/P4_casscf_bagel/P4_casscf_ang.json @@ -0,0 +1,27 @@ +{ "bagel" : [ + +{ + "title" : "molecule", + "basis" : "/home/ccsong/opt/bagel/current/share/cc-pvdz.json", + "df_basis" : "/home/ccsong/opt/bagel/current/share/cc-pvdz-jkfit.json", + "angstrom" : true, + "geometry" : [ + { "atom" : "P", "xyz" : [ 0.2238578355, 1.1309802898, 1.0989460181 ] }, + { "atom" : "P", "xyz" : [ 0.8151038290, -0.3588758134, 2.5788731692 ] }, + { "atom" : "P", "xyz" : [ -0.7631190474, -0.8141132093, 1.1430397092 ] }, + { "atom" : "P", "xyz" : [ 1.3075201574, -0.6565274109, 0.4744772210 ] } + ] +}, + +{ + "title" : "force", + "method" : [ { + "title" : "casscf", + "nact" : 12, + "nclosed" : 24, + "active" : [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36], + "maxiter" : 200, + "thresh" : 1e-10 + } ] +} +]} diff --git a/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang.log b/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang.log new file mode 100644 index 00000000..c6831c1a --- /dev/null +++ b/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang.log @@ -0,0 +1,122 @@ +geometric-optimize called with the following command line: +/home/leeping/opt/miniconda3/envs/main/bin/geometric-optimize --engine bagel P4_casscf_ang.json + + ())))))))))))))))/ + ())))))))))))))))))))))))), + *))))))))))))))))))))))))))))))))) + #, ()))))))))/ .)))))))))), + #%%%%, ()))))) .))))))))* + *%%%%%%, )) .. ,))))))). + *%%%%%%, ***************/. .))))))) + #%%/ (%%%%%%, /*********************. ))))))) + .%%%%%%# *%%%%%%, *******/, **********, .)))))) + .%%%%%%/ *%%%%%%, ** ******** .)))))) + ## .%%%%%%/ (%%%%%%, ,****** /))))) + %%%%%% .%%%%%%# *%%%%%%, ,/////. ****** )))))) + #% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,))))) + #%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))). + #%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))), + #%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))). + ## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,))))) + #%%%%# /%%%%%%/ (%%%%%% /)/)// ****** )))))) + ## .%%%%%%/ (%%%%%%, ******* )))))) + .%%%%%%/ *%%%%%%, **. /******* .)))))) + *%%%%%%/ (%%%%%% ********/*..,*/********* *)))))) + #%%/ (%%%%%%, *********************/ ))))))) + *%%%%%%, ,**************/ ,))))))/ + (%%%%%% () )))))))) + #%%%%, ()))))) ,)))))))), + #, ()))))))))) ,)))))))))). + ()))))))))))))))))))))))))))))))/ + ())))))))))))))))))))))))). + ())))))))))))))), + +-=#  geomeTRIC started. Version: 1.0.1+228.g3731120  #=- +Current date and time: 2024-06-06 11:47:14 + #========================================================# +#|  Arguments passed to driver run_optimizer():  |# +#========================================================# +engine bagel +input P4_casscf_ang.json +---------------------------------------------------------- +Bagel engine selected. Expecting Bagel input for gradient calculation. +Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii +12 internal coordinates being used (instead of 12 Cartesians) +Internal coordinate system (atoms numbered from 1): +Distance 1-2 +Distance 1-3 +Distance 1-4 +Distance 2-3 +Distance 2-4 +Distance 3-4 +Angle 2-1-3 +Angle 2-1-4 +Angle 3-1-4 +Angle 1-2-3 +Angle 1-2-4 +Angle 3-2-4 +Angle 1-3-2 +Angle 1-3-4 +Angle 2-3-4 +Angle 1-4-2 +Angle 1-4-3 +Angle 2-4-3 +Dihedral 3-1-2-4 +Dihedral 4-1-2-3 +Dihedral 2-1-3-4 +Dihedral 4-1-3-2 +Dihedral 2-1-4-3 +Dihedral 3-1-4-2 +Dihedral 1-2-3-4 +Dihedral 4-2-3-1 +Dihedral 1-2-4-3 +Dihedral 3-2-4-1 +Dihedral 1-3-4-2 +Dihedral 2-3-4-1 +Translation-X 1-4 +Translation-Y 1-4 +Translation-Z 1-4 +Rotation-A 1-4 +Rotation-B 1-4 +Rotation-C 1-4 + : 6 + : 12 + : 12 + : 1 + : 1 + : 1 + : 1 + : 1 + : 1 +> ===== Optimization Info: ==== +> Job type: Energy minimization +> Maximum number of optimization cycles: 300 +> Initial / maximum trust radius (Angstrom): 0.100 / 0.300 +> Convergence Criteria: +> Will converge when all 5 criteria are reached: +> |Delta-E| < 1.00e-06 +> RMS-Grad < 3.00e-04 +> Max-Grad < 4.50e-04 +> RMS-Disp < 1.20e-03 +> Max-Disp < 1.80e-03 +> === End Optimization Info === +Step 0 : Gradient = 4.614e-02/4.614e-02 (rms/max) Energy = -1363.1059640900 +Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.94892e-01 1.94902e-01 1.99904e-01 +Step 1 : Displace = 3.053e-02/3.054e-02 (rms/max) Trust = 1.000e-01 (=) Grad = 1.818e-02/1.819e-02 (rms/max) E (change) = -1363.1132801200 (-7.316e-03) Quality = 1.374 +Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.94887e-01 1.94892e-01 1.94902e-01 +Step 2 : Displace = 1.986e-02/1.986e-02 (rms/max) Trust = 1.414e-01 (+) Grad = 2.802e-03/2.803e-03 (rms/max) E (change) = -1363.1148298400 (-1.550e-03) Quality = 1.135 +Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.94887e-01 1.94892e-01 1.94902e-01 +Step 3 : Displace = 3.618e-03/3.619e-03 (rms/max) Trust = 2.000e-01 (+) Grad = 2.152e-04/2.154e-04 (rms/max) E (change) = -1363.1148709600 (-4.112e-05) Quality = 1.073 +Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.94887e-01 1.94892e-01 1.94902e-01 +Step 4 : Displace = 3.009e-04/3.013e-04 (rms/max) Trust = 2.828e-01 (+) Grad = 2.865e-06/2.999e-06 (rms/max) E (change) = -1363.1148712000 (-2.400e-07) Quality = 0.981 +Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.94887e-01 1.94892e-01 1.94902e-01 +Converged! =D + + #==========================================================================# + #| If this code has benefited your research, please support us by citing: |# + #| |# + #| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |# + #| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |# + #| http://dx.doi.org/10.1063/1.4952956 |# + #==========================================================================# + Time elapsed since start of run_optimizer: 1906.633 seconds diff --git a/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang_optim.xyz b/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang_optim.xyz new file mode 100644 index 00000000..7e256f03 --- /dev/null +++ b/examples/1-simple-examples/P4_casscf_bagel/saved/2024-06-07/P4_casscf_ang_optim.xyz @@ -0,0 +1,30 @@ +4 +Iteration 0 Energy -1363.10596409 +P 0.2238578355 1.1309802898 1.0989460181 +P 0.8151038290 -0.3588758134 2.5788731692 +P -0.7631190474 -0.8141132093 1.1430397092 +P 1.3075201574 -0.6565274109 0.4744772210 +4 +Iteration 1 Energy -1363.11328012 +P 0.2199265944 1.1608241814 1.0937984221 +P 0.8246799302 -0.3630969191 2.6075499860 +P -0.7896105087 -0.8287271310 1.1389210150 +P 1.3283667586 -0.6675362708 0.4550666944 +4 +Iteration 2 Energy -1363.11482984 +P 0.2173693822 1.1802371109 1.0904500778 +P 0.8309084669 -0.3658428839 2.6262026041 +P -0.8068423808 -0.8382330794 1.1362425736 +P 1.3419273131 -0.6746972871 0.4424408620 +4 +Iteration 3 Energy -1363.11487096 +P 0.2169034436 1.1837741460 1.0898398386 +P 0.8320427307 -0.3663436186 2.6296000583 +P -0.8099817908 -0.8399648697 1.1357554196 +P 1.3443984035 -0.6760017921 0.4401408011 +4 +Iteration 4 Energy -1363.11487120 +P 0.2168646647 1.1840685331 1.0897889677 +P 0.8321367852 -0.3663855293 2.6298822117 +P -0.8102429401 -0.8401088929 1.1357153503 +P 1.3446042773 -0.6761102499 0.4399495925