BAGEL engine

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
+  } ]
+}
+]}

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
+<class 'geometric.internal.Distance'> : 6
+<class 'geometric.internal.Angle'> : 12
+<class 'geometric.internal.Dihedral'> : 12
+<class 'geometric.internal.TranslationX'> : 1
+<class 'geometric.internal.TranslationY'> : 1
+<class 'geometric.internal.TranslationZ'> : 1
+<class 'geometric.internal.RotationA'> : 1
+<class 'geometric.internal.RotationB'> : 1
+<class 'geometric.internal.RotationC'> : 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

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

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:

geometric/errors.py

@@ -108,6 +108,9 @@ class GaussianEngineError(EngineError):
 class QCEngineAPIEngineError(EngineError):
     pass
 
+class BagelEngineError(EngineError):
+    pass
+
 class ConicalIntersectionEngineError(EngineError):
     pass
 

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')