ASH bugfixes

README.md

@@ -11,13 +11,16 @@ Interfaces to popular QM codes: ORCA, xTB, PySCF, MRCC, ccpy, Psi4, Dalton, CFou
 Excellent environment for writing simple or complex computational chemistry workflows.
 
 **Citation**
-If ASH is useful in your research please cite us: `ASH: a Multi-scale, Multi-theory Modeling program <https://onlinelibrary.wiley.com/doi/10.1002/jcc.70359>`_
+
+If ASH is useful in your research please cite us: 
+`ASH: a Multi-scale, Multi-theory Modeling program <https://onlinelibrary.wiley.com/doi/10.1002/jcc.70359>`_
 R. Bjornsson, J. Comput. Chem 2026, 47, e70359.
 
 **In case of problems:**
 Please open an issue on Github and we will try to fix any problems as soon as possible.
 
 **Installation:**
+
 See https://ash.readthedocs.io/en/latest/setup.html for detailed installation instructions.
 A proper ASH installation should usually be done in a conda/mamba environment together with the OpenMM library.
 

ash/__init__.py

@@ -80,7 +80,7 @@
 # Spinprojection
 from .modules.module_spinprojection import SpinProjectionTheory
 # HybridTheory: DualTheory and WrapTheory
-from .modules.module_hybridtheory import DualTheory,WrapTheory
+from .modules.module_hybridtheory import DualTheory,WrapTheory,FractTheory
 #ONIOM
 from .modules.module_oniom import ONIOMTheory
 

ash/functions/functions_optimization.py

@@ -561,7 +561,6 @@ def periodic_optimizer_alternating(fragment=None, theory=None, rate=0.5, maxiter
 
 # Cartesian-based periodic cell optimizer
 
-
 # Wrapper function around Cart_optimizer_class
 def Cart_optimizer(fragment=None, theory=None, rate=2.0, 
                                 scaling_rate_cell=1.0, maxiter=50, 
@@ -1084,6 +1083,9 @@ def compute_bfgs_step(self, current_grad, current_coords):
             # If curvature is bad, reset the Hessian to Identity to avoid exploding
             print("BFGS: Curvature condition not met, resetting Hessian.")
             self.Hess_inv = np.eye(n) * self.rate
+            # Also reset history to avoid stale s/y being used next step
+            self.g_old = g
+            self.x_old = x
 
         # 4. COMPUTE STEP
         # p = -Hess_inv * g
@@ -1164,6 +1166,7 @@ def compute_step(self,gradient,currcoords):
             effective_gradient = gradient * rate_mask
         else:
             effective_gradient = gradient
+
         # Calculate delta step (in Bohrs)
         if self.step_algo.lower() =="sd":
             print("Taking steepest descent step")
@@ -1186,7 +1189,7 @@ def compute_step(self,gradient,currcoords):
             delta_au = nesterov_update
         elif self.step_algo.lower() == "bfgs":
             print("Taking BFGS step")
-            delta_au = self.compute_bfgs_step(gradient, currcoords)
+            delta_au = self.compute_bfgs_step(effective_gradient, currcoords)
         elif self.step_algo.lower() == "cg":
             print("Taking conjugate gradient step")
             if self.iteration == 0:
@@ -1205,6 +1208,11 @@ def compute_step(self,gradient,currcoords):
             print("Unknown step_algo")
             ashexit()
 
+        # If cell is frozen, explicitly zero those rows in the step too,
+        # so accumulated BFGS state can never leak cell movement back in
+        if self.PBC and self.scaling_rate_cell == 0.0:
+            delta_au[-3:] = 0.0
+
         return delta_au
 
     def run(self, theory=None, fragment=None, constraints=None, charge=None, mult=None):
@@ -1463,4 +1471,4 @@ def run(self, theory=None, fragment=None, constraints=None, charge=None, mult=No
             #currcoords = currcoords_new_ang
 
         if iteration == self.maxiter-1:
-            print("Number of max iterations reached without reaching convergence. Sad...")
+            print("Number of max iterations reached without reaching convergence. Sad...")# Also reset history to avoid stale s/y being used next step

ash/interfaces/interface_CP2K.py

@@ -71,7 +71,7 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel
                 print("potential_dict keyword is required")
                 ashexit()
             if functional is None:
-                print("functional keyword is required for PW andd GPW ")
+                print("functional keyword is required for PW and GPW ")
                 ashexit()
         else:
             print("This is a CP2K xTB theory")
@@ -448,7 +448,7 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el
                              basis_dict=self.basis_dict, potential_dict=self.potential_dict,
                              basis_method=self.basis_method, wavelet_scf_type=self.wavelet_scf_type,
                              functional=self.functional, restartfile=None, mgrid_commensurate=True,
-                             Grad=Grad, filename='cp2k', charge=charge, mult=mult,
+                             Grad=Grad, filename=self.filename, charge=charge, mult=mult,
                              coordfile=system_xyzfile,
                              stress_tensor=self.stress_tensor, stress_tensor_algo=self.stress_tensor_algo,
                              user_input_dft=self.user_input_dft, vdwpotential=self.vdwpotential,
@@ -497,7 +497,7 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el
                              basis_dict=self.basis_dict, potential_dict=self.potential_dict,
                              basis_method=self.basis_method, wavelet_scf_type=self.wavelet_scf_type,
                              functional=self.functional, restartfile=None,
-                             Grad=Grad, filename='cp2k', charge=charge, mult=mult,
+                             Grad=Grad, filename=self.filename, charge=charge, mult=mult,
                              stress_tensor=self.stress_tensor, stress_tensor_algo=self.stress_tensor_algo,
                              user_input_dft=self.user_input_dft, vdwpotential=self.vdwpotential,
                              kpoint_settings=self.kpoint_settings,
@@ -1068,4 +1068,57 @@ def stress_to_cell_gradient(lattice_matrix, stress_tensor):
     # compute gradient
     grad = -1*V * sigma @ np.linalg.inv(h).T
 
-    return grad
+    return grad
+
+
+def grabatomcharges_CP2K(outputfile, chargemodel="Mulliken"):
+
+    with open(outputfile, "r") as f:
+        lines = f.readlines()
+
+        chargemodel = chargemodel.lower()
+        charge_sections = {"mulliken": "Mulliken Population Analysis", "hirshfeld": "Hirshfeld Charges"}
+
+    if chargemodel not in charge_sections:
+        print(f"Unknown CP2K charge model: {chargemodel}")
+        ashexit()
+
+    section = charge_sections[chargemodel]
+
+    # Find last occurrence, since output may contain multiple calculations # Check with multiple single points
+    start = None
+    for i, line in enumerate(lines):
+        if section.lower() in line.lower():
+            start = i
+
+    if start is None:
+        print(f"Could not find {section} in CP2K output")
+        ashexit()
+
+    charges = []
+
+    for line in lines[start + 1:]:
+        parts = line.split()
+
+        if not parts:
+            continue
+
+        # Atom number must be first column
+        try:
+            int(parts[0])
+        except ValueError:
+            if charges:
+                break
+            continue
+
+        # Last column = Net charge
+        try:
+            charge = float(parts[-1])
+        except ValueError:
+            continue
+
+        charges.append(charge)
+
+    print(f"Grabbed {len(charges)} {chargemodel} charges from CP2K")
+
+    return charges

ash/interfaces/interface_ORCA.py

@@ -372,12 +372,34 @@ def get_dipole_moment(self):
         dm = grab_dipole_moment(self.filename+'.out')
         print("Dipole moment:", dm)
         return dm
+
     def get_polarizability_tensor(self):
         print("here")
         print("self.filename+'.out':", self.filename+'.out')
         polarizability,diag_pz = grab_polarizability_tensor(self.filename+'.out')
         print("polarizability:", polarizability)
         return polarizability
+
+    # Get density matrix. Assuming dm has been defined
+    def get_density_matrix(self, dmoption=None):
+        print("Calling get_density_matrix")
+        # Create JSON-file from ORCA-GBW and density-files
+        jsonfile = create_ORCA_json_file(self.filename+'.gbw', format="json", fock_matrix=True)
+        #Read the JSON-file
+        data = read_ORCA_json_file(jsonfile)
+        fock = np.array(data.get("F-Matrix", None)[0])
+        #Get densities from data dictionary (from read_ORCA_json_file)
+        DMs = get_densities_from_ORCA_json(data)
+        print("DMs:", DMs)
+        if dmoption is None:
+            print("No dmoption provided to get_density_matrix. Defaulting to 'scfp' density")
+            dmoption = "scfp"
+        #Grab ORCA density from jsonfile or data-dictionary. Returns DM_AO,C,S, MO_occs, MO_energies, AO_basis, AO_order
+        DM_AO,C,S, MO_occs, MO_energies, AO_basis, AO_order = grab_ORCA_wfn(jsonfile=jsonfile, 
+                                                                            density=dmoption)
+        print("DM_AO:", DM_AO)
+        return DM_AO, fock
+
     # Run function. Takes coords, elems etc. arguments and computes E or E+G.
     def run(self, current_coords=None, charge=None, mult=None, current_MM_coords=None, MMcharges=None, qm_elems=None, mm_elems=None,
             elems=None, Grad=False, Hessian=False, PC=False, numcores=None, label=None):
@@ -3059,7 +3081,8 @@ def create_GBW_from_json_file(jsonfile, orcadir=None):
 #Using orca_2json to create JSON file from ORCA GBW file
 #Format options: json, bson, ubjson, msgpack
 def create_ORCA_json_file(file, orcadir=None, format="json", basis_set=True, mo_coeffs=True, one_el_integrals=True,
-                          two_el_integrals=False, two_el_integrals_type="ALL", dipole_integrals=False, full_int_transform=False):
+                          two_el_integrals=False, two_el_integrals_type="ALL", dipole_integrals=False, 
+                          full_int_transform=False, fock_matrix=False):
     print("create_ORCA_json_file")
     orcadir = check_ORCA_location(orcadir, modulename="create_ORCA_json_file")
     #orcafile_basename = file.split('.')[0]
@@ -3071,6 +3094,7 @@ def create_ORCA_json_file(file, orcadir=None, format="json", basis_set=True, mo_
     two_el_integrals_line=""
     basis_set_line=""
     mo_coeff_line=""
+    fock_matrix_line=""
     #NOTE: problems with FullTrafo (orca_2json crashes)
     if full_int_transform is True:
         full_transform_integrals_line="\"FullTrafo\": true,"
@@ -3093,6 +3117,9 @@ def create_ORCA_json_file(file, orcadir=None, format="json", basis_set=True, mo_
             two_el_integrals_line=f"\"2elIntegrals\": [\"MO_PQRS\", \"MO_PRQS\"],"
         else:
             two_el_integrals_line=f"\"2elIntegrals\": [\"MO_{two_el_integrals_type}\"],"
+    if fock_matrix is True:
+        print("Requesting printout of Fock matrix")
+        fock_matrix_line="\"FockMatrix\": [\"F\"],"
     if mo_coeffs is True:
         print("Requesting printout of MO coefficients")
         mo_coeff_line="\"MOCoefficients\": true,"
@@ -3105,6 +3132,7 @@ def create_ORCA_json_file(file, orcadir=None, format="json", basis_set=True, mo_
 {prop_1e_integrals_line}
 {two_el_integrals_line}
 {full_transform_integrals_line}
+{fock_matrix_line}
 "Densities": ["all"],
 "JSONFormats": ["{format}"]
 }}

ash/interfaces/interface_OpenMM.py

@@ -3952,51 +3952,60 @@ def __init__(self, fragment=None, theory=None, charge=None, mult=None, timestep=
     def set_sim_reporters(self,simulation,restart=False):
         import openmm
         #CheckpointReporter
-        print("Creating CheckpointReporter that will write a restartable checkpointfile every X steps")
-        checkpointfilename='OpenMM_MD.chk'
-        simulation.reporters.append(openmm.app.CheckpointReporter(checkpointfilename, self.traj_frequency*1))
+        if restart is False:
+            print("Creating CheckpointReporter that will write a restartable checkpointfile every X steps")
+            checkpointfilename='OpenMM_MD.chk'
+            simulation.reporters.append(openmm.app.CheckpointReporter(checkpointfilename, self.traj_frequency*1))
         #StateDataReporter
-        print("Creating StateDataReporter that will write to stdout")
-        statedatareporter_stdout=openmm.app.StateDataReporter(stdout, self.traj_frequency, step=True, time=True,
-                                                    potentialEnergy=True, kineticEnergy=True, volume=self.volume,
-                                                    density=self.density, temperature=True, separator=',')
-        simulation.reporters.append(statedatareporter_stdout)
+        if restart is False:
+            print("Creating StateDataReporter that will write to stdout")
+            statedatareporter_stdout=openmm.app.StateDataReporter(stdout, self.traj_frequency, step=True, time=True,
+                                                        potentialEnergy=True, kineticEnergy=True, volume=self.volume,
+                                                        density=self.density, temperature=True, separator=',')
+
+            simulation.reporters.append(statedatareporter_stdout)
         #Another reporter for writing to file
         if self.dataoutputoption != stdout:
-            print("Creating StateDataReporter that will write to file:", self.datafilename)
-            print("restart:", restart)
-            statedatareporter_file=openmm.app.StateDataReporter(self.dataoutputoption, self.traj_frequency, step=True, time=True,
-                                                            potentialEnergy=True, kineticEnergy=True, volume=self.volume,
-                                                            density=self.density, temperature=True, separator=',', append=restart)
-            simulation.reporters.append(statedatareporter_file)
-            self.dataoutputoption = open(self.datafilename,'a')
+
+            if restart is False:
+                print("Creating StateDataReporter that will write to file:", self.datafilename)
+                print("restart:", restart)
+                statedatareporter_file=openmm.app.StateDataReporter(self.dataoutputoption, self.traj_frequency, step=True, time=True,
+                                                                potentialEnergy=True, kineticEnergy=True, volume=self.volume,
+                                                                density=self.density, temperature=True, separator=',', append=restart)
+                simulation.reporters.append(statedatareporter_file)
+                self.dataoutputoption = open(self.datafilename,'a')
 
         #TODO: See if this can be made to work for simulations with step-by-step
         if self.trajectory_file_option == 'PDB':
-            simulation.reporters.append(
-                openmm.app.PDBReporter(self.trajfilename+'.pdb', self.traj_frequency,
-                                                         enforcePeriodicBox=self.enforcePeriodicBox))
+            if restart is False:
+                simulation.reporters.append(
+                    openmm.app.PDBReporter(self.trajfilename+'.pdb', self.traj_frequency,
+                                                            enforcePeriodicBox=self.enforcePeriodicBox))
         elif self.trajectory_file_option == 'DCD':
             #Note: using append keyword here if restarting
             # Check first if file exists for restart (OpenMM errors otherwise)
             if restart is True and os.path.isfile(f"{self.trajfilename}.dcd") is False:
                 print("Warning: restart option was active but trajectory file not existing. Will create new file")
                 restart=False
 
-            simulation.reporters.append(
-                openmm.app.DCDReporter(self.trajfilename+'.dcd', self.traj_frequency, append=restart,
-                                                         enforcePeriodicBox=self.enforcePeriodicBox))
-            print('DCDReporter added')
+            if restart is False:
+                simulation.reporters.append(
+                    openmm.app.DCDReporter(self.trajfilename+'.dcd', self.traj_frequency, append=restart,
+                                                            enforcePeriodicBox=self.enforcePeriodicBox))
+                print('DCDReporter added')
         elif self.trajectory_file_option == 'NetCDFReporter':
             print("NetCDFReporter traj format selected. This requires mdtraj. Importing.")
             mdtraj = MDtraj_import()
-            simulation.reporters.append(
-                mdtraj.reporters.NetCDFReporter(self.trajfilename+'.nc', self.traj_frequency))
+            if restart is False:
+                simulation.reporters.append(
+                    mdtraj.reporters.NetCDFReporter(self.trajfilename+'.nc', self.traj_frequency))
         elif self.trajectory_file_option == 'HDF5Reporter':
             print("HDF5Reporter traj format selected. This requires mdtraj. Importing.")
             mdtraj = MDtraj_import()
-            simulation.reporters.append(
-                mdtraj.reporters.HDF5Reporter(self.trajfilename+'.lh5', self.traj_frequency,
+            if restart is False:
+                simulation.reporters.append(
+                    mdtraj.reporters.HDF5Reporter(self.trajfilename+'.lh5', self.traj_frequency,
                                               enforcePeriodicBox=self.enforcePeriodicBox))
         elif self.trajectory_file_option =="XYZ":
             print("XYZ trajectory format selected (not available for classical MD). Warning: not very fast")
@@ -4010,7 +4019,8 @@ def set_sim_reporters(self,simulation,restart=False):
         if self.force_file_option != None:
             print("ForceReporter traj format selected.")
             #exit()
-            simulation.reporters.append(ForceReporter(self.trajfilename + '_force.txt', self.traj_frequency, atomic_units=self.atomic_units_force_reporter))
+            if restart is False:
+                simulation.reporters.append(ForceReporter(self.trajfilename + '_force.txt', self.traj_frequency, atomic_units=self.atomic_units_force_reporter))
         if self.energy_file_option != None:
             print("Energyfile  selected.")
             try: