diff --git a/README.md b/README.md index 7b4eb651a..4bad11fc2 100644 --- a/README.md +++ b/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 `_ + +If ASH is useful in your research please cite us: +`ASH: a Multi-scale, Multi-theory Modeling program `_ 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. diff --git a/ash/__init__.py b/ash/__init__.py index ac5abb93a..9b82be0ac 100644 --- a/ash/__init__.py +++ b/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 diff --git a/ash/functions/functions_optimization.py b/ash/functions/functions_optimization.py index c43c6e353..47bcbe223 100644 --- a/ash/functions/functions_optimization.py +++ b/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...") \ No newline at end of file + print("Number of max iterations reached without reaching convergence. Sad...")# Also reset history to avoid stale s/y being used next step diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 682ccc5f1..399edaef7 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/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 \ No newline at end of file + 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 \ No newline at end of file diff --git a/ash/interfaces/interface_ORCA.py b/ash/interfaces/interface_ORCA.py index 312b2b58d..8181a08e4 100644 --- a/ash/interfaces/interface_ORCA.py +++ b/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}"] }} diff --git a/ash/interfaces/interface_OpenMM.py b/ash/interfaces/interface_OpenMM.py index f851fa7c7..e802ad07b 100644 --- a/ash/interfaces/interface_OpenMM.py +++ b/ash/interfaces/interface_OpenMM.py @@ -3952,30 +3952,36 @@ 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) @@ -3983,20 +3989,23 @@ def set_sim_reporters(self,simulation,restart=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: diff --git a/ash/interfaces/interface_block.py b/ash/interfaces/interface_block.py index 93e52fef6..7b735cca0 100644 --- a/ash/interfaces/interface_block.py +++ b/ash/interfaces/interface_block.py @@ -32,7 +32,7 @@ def __init__(self, blockdir=None, pyscftheoryobject=None, blockversion='Block2', block_parallelization='OpenMP', numcores=1, hybrid_num_mpi_procs=None, hybrid_num_threads=None, FIC_MRCI=False, SC_NEVPT2_Wick=False, IC_NEVPT2=False, DMRG_DoRDM=False, DMRG_DoRDM2=False, SC_NEVPT2=False, SC_NEVPT2_Mcompression=None, label="Block", print_WF_coeffs=False, - groupname=None, orbsym=None): + groupname=None, orbsym=None, restart=False): self.theorynamelabel="Block" self.theorytype="QM" @@ -106,6 +106,8 @@ def __init__(self, blockdir=None, pyscftheoryobject=None, blockversion='Block2', ashexit() self.orbsym=orbsym + # RESTART + self.restart=restart # untested #SETTING NUMCORES by setting prefix self.block_parallelization=block_parallelization @@ -451,19 +453,15 @@ def setup_DMRG_job(self,verbose=5, rdmoption=None): if self.macroiter == 0: print("This is single-iteration CAS-CI via pyscf and DMRG") #Creating pyscf CAS-CI object and setting fcisolver to DMRGCI - #print("self.pyscftheoryobject.mol:", self.pyscftheoryobject.mol) - #print("self.pyscftheoryobject.mol:", self.pyscftheoryobject.mol.__dict__) - #print("----------------") - #print("self.pyscftheoryobject.mf.mol:", self.pyscftheoryobject.mf.mol) - #print("self.pyscftheoryobject.mf.mol:", self.pyscftheoryobject.mf.mol.__dict__) - #exit() - #print("self.pyscftheoryobject.mf", self.pyscftheoryobject.mf) self.mch = self.pyscf.mcscf.CASCI(self.pyscftheoryobject.mf, self.norb, self.nelec) #print("self.mch:", self.mch) - self.mch.fcisolver = self.dmrgscf.DMRGCI(self.pyscftheoryobject.mol, maxM=self.maxM, tol=self.tol) - #print("self.mch.fcisolver:", self.mch.fcisolver) - #print("self.mch.fcisolver wfnsym:", self.mch.fcisolver.wfnsym) - #print("1self.mch.fcisolver.groupname:", self.mch.fcisolver.groupname) + + self.mch.fcisolver = self.dmrgscf.DMRGCI(self.pyscftheoryobject.mol, + maxM=self.maxM, tol=self.tol) + if self.restart is True: + print("Restarting DMRG-CASCI job.") + self.mch.fcisolver.restart = True + if self.groupname is not None: print("Setting groupname for mch.fcisolver and orbsym in mch") self.mch.fcisolver.groupname=self.groupname @@ -473,21 +471,15 @@ def setup_DMRG_job(self,verbose=5, rdmoption=None): if self.mch.mol.groupname == "N/A": self.mch.mol.groupname = 'C1' - #print("2self.mch.fcisolver.groupname:", self.mch.fcisolver.groupname) - #print("self.mch.orbsym:", self.mch.orbsym) - #print("self.mch.mol:", self.mch.mol) - #print("self.mch.mol dict:", self.mch.mol.__dict__) - #print("self.mch.mol.groupname:", self.mch.mol.groupname) - #self.mch = self.pyscf.mcscf.CASCI(self.pyscftheoryobject.mf,self.norb, self.nelec) - #self.mch = self.dmrgscf.DMRGCI(self.pyscftheoryobject.mf,self.norb, self.nelec, maxM=self.maxM, tol=self.tol) - #self.mch = self.dmrgscf.DMRGSCF(self.pyscftheoryobject.mf, self.norb, self.nelec, maxM=self.maxM, tol=self.tol) - #print("Turning off canonicalization step in mcscf object") - #self.mch.canonicalization = False - #self.mch.natorb = True else: print("This is CASSCF via pyscf and DMRG (orbital optimization)") # - self.mch = self.dmrgscf.DMRGSCF(self.pyscftheoryobject.mf,self.norb, self.nelec, maxM=self.maxM, tol=self.tol) + + self.mch = self.dmrgscf.DMRGSCF(self.pyscftheoryobject.mf,self.norb, self.nelec, + maxM=self.maxM, tol=self.tol, restart=self.restart) + if self.restart is True: + print("Restarting DMRG-CASSCF job.") + self.mch.fcisolver.restart = True self.mch.canonicalization = True self.mch.natorb = True self.mch.chkfile = f"DMRG-CASSCF_{self.maxM}.chk" @@ -623,12 +615,16 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el qm_elems = elems # Cleanup before run. - self.cleanup() - - # Run PySCF to get integrals and MOs. This would probably only be an SCF - if self.Block_direct != True: - print("Now running input PySCFTheory object") - self.pyscftheoryobject.run(current_coords=current_coords, elems=qm_elems, charge=charge, mult=mult) + if self.restart: + print("Restarting previous Block job. Will not cleanup") + else: + print("This is not a restarted job. Will cleanup.") + self.cleanup() + # Run PySCF to get integrals and MOs. This would probably only be an SCF + if self.Block_direct != True: + print("Block_direct is not True. Will run PySCF to get integrals and MOs") + print("Now running input PySCFTheory object") + self.pyscftheoryobject.run(current_coords=current_coords, elems=qm_elems, charge=charge, mult=mult) # Get frozen-core if self.frozencore is True: @@ -647,13 +643,14 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el print("Running Block directly") self.write_inputfile(mult) self.call_block_directly() + elif self.restart is True: + print("Restarting previous Block job. Will not setup initial orbitals") else: if self.runcalls == 1: print("First runcall.") print("Doing initial orbital setup") # Natural orbital if self.initial_orbitals != "HF": - print mo_coeffs, occupations = self.setup_initial_orbitals(elems) #Returns mo-coeffs and occupations of initial orbitals print("Doing active space setup") self.setup_active_space(occupations=occupations) #This will define self.norb and self.nelec active space diff --git a/ash/interfaces/interface_dlfind.py b/ash/interfaces/interface_dlfind.py index caed45192..da619f526 100644 --- a/ash/interfaces/interface_dlfind.py +++ b/ash/interfaces/interface_dlfind.py @@ -15,7 +15,7 @@ from ash.modules.module_coords_PBC import write_CIF_file, write_XSF_file, write_POSCAR_file, cell_vectors_to_params, cell_volume, align_to_standard_orientation from ash.modules.module_theory import NumGradclass from ash.modules.module_results import ASH_Results -from ash.modules.module_freq import NumFreq,AnFreq,calc_hessian_xtb +from ash.modules.module_freq import NumFreq,AnFreq, approximate_full_Hessian_from_smaller,calc_hessian_xtb from ash.modules.module_QMMM import QMMMTheory from ash.modules.module_oniom import ONIOMTheory @@ -335,7 +335,14 @@ def hess_func(coords): hessatoms=self.numfreq_hessatoms,force_projection=self.numfreq_force_projection, runmode='serial', numcores=self.theory.numcores) - hessian = result_freq.hessian + + if self.numfreq_hessatoms is not None and len(self.numfreq_hessatoms) < self.fragment.numatoms: + print("A partial Hessian was computed. Approximating full Hessian from smaller Hessian") + hessian = approximate_full_Hessian_from_smaller(fragment,result_freq.hessian, + self.numfreq_hessatoms, large_atomindices=self.fragment.allatoms, + restHessian=None) + else: + hessian = result_freq.hessian elif self.hessian_choice == "anfreq": print("AnFreq option requested") result_freq = AnFreq(theory=self.theory, fragment=self.fragment, printlevel=0) diff --git a/ash/interfaces/interface_pyscf.py b/ash/interfaces/interface_pyscf.py index f065a4c6c..fc367adf9 100644 --- a/ash/interfaces/interface_pyscf.py +++ b/ash/interfaces/interface_pyscf.py @@ -877,6 +877,17 @@ def run_population_analysis(self, mf, unrestricted=False, dm=None, type='Mullike print(f"{label} Mulliken spin pops:", mulliken_spinpopulations[1]) return + # Get density matrix. Assuming dm has been defined + def get_density_matrix(self): + + # Check if self.dm is defined + if not hasattr(self, 'dm'): + print("Density matrix not defined. Calculating from mf object") + dm = self.mf.make_rdm1() + else: + dm = self.dm + return dm + def run_stability_analysis(self): def stableprint(stable_i,stable_e): if stable_i is True: diff --git a/ash/interfaces/interface_xtb.py b/ash/interfaces/interface_xtb.py index 09d254b4a..e8b54d3bd 100644 --- a/ash/interfaces/interface_xtb.py +++ b/ash/interfaces/interface_xtb.py @@ -547,6 +547,12 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el # Check if finished. Grab energy, gradient, pcgradient, cellgradient if Grad is True: self.energy,self.grad=xtbgradientgrab(num_qmatoms) + + # Fix for possible problem getting energy from gradient file (can happen for big systems) + if self.energy is None: + # Grabbing energy from main outputfile instead + self.energy=xtbfinalenergygrab(self.filename+'.out') + if self.periodic: self.cell_gradient = grab_latticegrad() print("cell_gradient:", self.cell_gradient) @@ -754,7 +760,12 @@ def xtbgradientgrab(numatoms): with open('gradient') as f: for line in reverse_lines(f): if ' cycle =' in line: - energy = float(line.split("SCF energy =")[1].split()[0]) + # Checking if energy is a number or asterisk (can happen for big systems) + en_entry = line.split("SCF energy =")[1].split()[0] + if '*' in en_entry: + energy = None + else: + energy = float(line.split("SCF energy =")[1].split()[0]) return energy, gradient if count==numatoms: grab=False diff --git a/ash/modules/module_freq.py b/ash/modules/module_freq.py index 9c0e512a4..5c621fe51 100644 --- a/ash/modules/module_freq.py +++ b/ash/modules/module_freq.py @@ -521,9 +521,9 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa # Use input masses if given, otherwise take from frament if hessatoms_masses is None: - print("allatoms:", allatoms) - print("hessatoms:", hessatoms) - print("fragment.list_of_masses:", fragment.list_of_masses) + #print("allatoms:", allatoms) + #print("hessatoms:", hessatoms) + #print("fragment.list_of_masses:", fragment.list_of_masses) hessmasses = ash.modules.module_coords.get_partial_list(allatoms, hessatoms, fragment.list_of_masses) else: hessmasses=hessatoms_masses @@ -1376,7 +1376,7 @@ def calc_model_Hessian_ORCA(fragment,model='Almloef'): # capping_atom_hessian_indices=[] #NOTE: Trans+rot projection off right now def approximate_full_Hessian_from_smaller(fragment, hessian_small, small_atomindices, large_atomindices=None, restHessian='zero', projection=False, - charge=None, mult=None, xtbmethod="GFN1"): + charge=None, mult=None, xtbmethod="GFN1", diagonalize_fullhessian=False): print("approximate_full_Hessian_from_smaller") print() write_hessian(hessian_small,hessfile="smallhessian") @@ -1426,7 +1426,7 @@ def approximate_full_Hessian_from_smaller(fragment, hessian_small, small_atomind fullhessian=np.zeros((hess_size,hess_size)) print("Initial fullhessian:", fullhessian) print("Number of Hessian elements:", fullhessian.size) - write_hessian(fullhessian,hessfile="initialfullhessian") + #write_hessian(fullhessian,hessfile="initialfullhessian") #Making sure hessian_small is np array hessian_small = np.array(hessian_small) @@ -1467,7 +1467,7 @@ def approximate_full_Hessian_from_smaller(fragment, hessian_small, small_atomind print("RestHessian is zero.") print("Intermediate fullhessian:", fullhessian) print("Size:", fullhessian.size) - write_hessian(fullhessian,hessfile="intermedfullhessian") + #write_hessian(fullhessian,hessfile="intermedfullhessian") #Large Hessian indices athessindices = [3*i+j for i in correct_small_atomindices for j in [0,1,2]] #Looping over and assigning small Hessian values to large @@ -1475,7 +1475,7 @@ def approximate_full_Hessian_from_smaller(fragment, hessian_small, small_atomind for s_j, j in enumerate(athessindices): fullhessian[i,j] = hessian_small[s_i,s_j] print("Final fullhessian:", fullhessian) - write_hessian(fullhessian,hessfile="intermedfullhessian_after_small_update") + #write_hessian(fullhessian,hessfile="intermedfullhessian_after_small_update") #NOTE: Diagonalizing full Hessian just to see #Checking for linearity. Determines how many Trans+Rot modes if detect_linear(coords=fragment.coords,elems=fragment.elems) is True: @@ -1483,11 +1483,14 @@ def approximate_full_Hessian_from_smaller(fragment, hessian_small, small_atomind else: TRmodenum=6 - print("Now diagonalizing full Hessian") - frequencies, normal_modes, evectors, mode_order = diagonalizeHessian(fragment.coords,fullhessian,usedfragment.masses,usedfragment.elems,TRmodenum=TRmodenum,projection=projection) - print("Size:", fullhessian.size) - print("Frequencies of full Hessian:", frequencies) - write_hessian(fullhessian,hessfile="Finalfullhessian") + if diagonalize_fullhessian is True: + print("Now diagonalizing full Hessian") + frequencies, normal_modes, evectors, mode_order = diagonalizeHessian(fragment.coords,fullhessian, + usedfragment.masses,usedfragment.elems,TRmodenum=TRmodenum,projection=projection) + print("Frequencies of full Hessian:", frequencies) + print("Size of full Hessian:", fullhessian.size) + + #write_hessian(fullhessian,hessfile="Finalfullhessian") return fullhessian diff --git a/ash/modules/module_hybridtheory.py b/ash/modules/module_hybridtheory.py index 1b4299b8c..12a6f7450 100644 --- a/ash/modules/module_hybridtheory.py +++ b/ash/modules/module_hybridtheory.py @@ -4,6 +4,8 @@ """ import time +from ash.interfaces.interface_pyscf import PySCFTheory +from ash.interfaces.interface_ORCA import ORCATheory import numpy as np import os import copy @@ -581,3 +583,138 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el return self.energy, self.gradient else: return self.energy + + +######################### +# FractTheory class +######################### + +class FractTheory: + def __init__(self,theory=None, eta=None, chargeA=None, chargeB=None, + multA=None, multB=None, DMinterpol=False, + resultfile="fract_theory_gaps.txt"): + + self.eta=eta + self.chargeA=chargeA + self.chargeB=chargeB + self.multA=multA + self.multB=multB + self.theory=theory + + self.theorynamelabel="Fract" + self.theorytype="QM" + self.analytic_hessian=False + self.printlevel=2 + self.DMinterpol=DMinterpol # density matrix interpolation or not + self.resultfile=resultfile + + # Check for compatibility of theory + + # pyscf + if isinstance(self.theory,PySCFTheory): + if self.DMinterpol: + # Defining theory that has density matrix + self.dmtheory=self.theory + # QM/MM + elif self.theory.theorytype == "QM/MM": + if self.DMinterpol: + print("Theory type is QM/MM and DMinterpol is True. Checking if QM theory is compatible") + if not isinstance(self.theory.qm_theory, PySCFTheory): + print("Error: QM-theory inside QMMMTheory is not yet compatible with FractTheory.") + ashexit() + self.dmtheory=self.qm_theory + # Wraptheory + elif self.theory.theorynamelabel == "WrapTheory": + if self.DMinterpol: + print("Theory is WrapTheory and DMinterpol is True. Checking if it contains a compatible QM-theory") + # Checking instance of all theories via list comprehension + ll = [isinstance(t, PySCFTheory) for t in self.theory.theories] + # Check if any of the theories are PySCFTheory + if not any(ll): + print("Error: No compatible QM-theory found inside WrapTheory.") + ashexit() + # Defining dm theory the first theory that is compatible + self.dmtheory=self.theory.theories[ll.index(True)] + + # Results (gaps) writing to file, here header. + with open(self.resultfile, 'w') as f: + f.write(f"# eta deltaE (Eh)\n") + + + def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None, mm_elems=None, + elems=None, Grad=False, PC=False, numcores=None, pe=False, potfile=None, restart=False, label=None, + charge=None, mult=None, Hessian=False): + + print("Running FractTheory") + + # Init + print("Now running inital state SCF") + res_init = self.theory.run(current_coords=current_coords, elems=elems, + charge=self.chargeA, mult=self.multA, Grad=Grad) + if Grad: + E_init = res_init[0] + G_init = res_init[1] + else: + E_init = res_init + if self.DMinterpol: + dm_init = self.dmtheory.get_density_matrix() + + # Ion + print("Now running final-state SCF") + res_final = self.theory.run(current_coords=current_coords, elems=elems, charge=self.chargeB, mult=self.multB, Grad=Grad) + if Grad: + E_fin = res_final[0] + G_fin = res_final[1] + else: + E_fin = res_final + if self.DMinterpol: + dm_Nm1 = self.dmtheory.get_density_matrix() + + # Printing E(init) and E(final) and gap + print("E(init)", E_init) + print("E(final)", E_fin) + deltaE = E_fin - E_init + print("deltaE:", deltaE) + + # Interpolated Energy + E_eta = (1 - self.eta) * E_init + self.eta * E_fin + print("E_eta:", E_eta) + + self.energy=E_eta + + # Interpolated Gradient + if Grad: + G_eta = (1 - self.eta) * G_init + self.eta * G_fin + print("G_eta:", G_eta) + self.gradient=G_eta + + # Write gap to file + with open(self.resultfile, 'a') as f: + f.write(f"{self.eta} {deltaE}\n") + + # Derive Fractional energy via interpolated DM + if self.DMinterpol is True: + # Derive interpolated density matrix + print("Deriving fractional density matrix using eta:", self.eta) + # Fractional electron ensemble density + dm_fne = (1 - self.eta) * dm_init + self.eta * dm_Nm1 + self.dm=dm_fne + print("Deriving new energy from energy of interpolated DM") + if isinstance(self.dmtheory,PySCFTheory): + interpol_DM_E = float(self.dmtheory.mf.energy_tot(dm_fne)) + print("Interpol-DM energy:", interpol_DM_E) + elif isinstance(self.dmtheory,ORCATheory): + ashexit() + # Derive new energy from FNE DM + # Todo: need Fock matrices, density matrices and 1-el integrals + else: + print("Error: Incompatible input theory for FractTheory. Currently supporting: PySCFTheory") + ashexit() + self.energy=interpol_DM_E + + print("FractTheory energy:", self.energy) + + if Grad is True: + return self.energy,self.gradient + else: + return self.energy \ No newline at end of file diff --git a/ash/modules/module_oniom.py b/ash/modules/module_oniom.py index a82724166..92f50631e 100644 --- a/ash/modules/module_oniom.py +++ b/ash/modules/module_oniom.py @@ -22,7 +22,8 @@ def __init__(self, theories_N=None, regions_N=None, regions_chargemult=None, fullregion_charge=None, fullregion_mult=None, fragment=None, label=None, chargeboundary_method="chargeshift", excludeboundaryatomlist=None, linkatom_method='ratio', linkatom_simple_distance=None, linkatom_forceproj_method="adv", - linkatom_ratio=0.723, linkatom_type='H', printlevel=2, numcores=1): + linkatom_ratio=0.723, linkatom_type='H', printlevel=2, numcores=1, + LL_full_theory=None): super().__init__() self.theorytype="ONIOM" self.printlevel=printlevel @@ -61,6 +62,17 @@ def __init__(self, theories_N=None, regions_N=None, regions_chargemult=None, self.fragment=fragment self.allatoms = self.fragment.allatoms self.theories_N=theories_N + # Optional separate Theory for the LL-on-Full-system calculation. + # If not provided, behavior is unchanged: theories_N[-1] is used for both + # the Full-region LL run AND the Model-region LL run (subtraction term). + # If provided, theories_N[-1] is used only as LL_model (region subtraction term), + # while LL_full_theory is used only for the Full-region run. + # This allows LL_full and LL_model to be different theory objects/codes. + self.LL_full_theory = LL_full_theory + if self.LL_full_theory is not None: + print(f"Separate LL_full_theory provided: {self.LL_full_theory.theorynamelabel}") + print("This will be used ONLY for the Full-region low-level calculation.") + print(f"theories_N[-1] ({theories_N[-1].theorynamelabel}) will still be used as LL_model (region subtraction term).") self.regions_N=regions_N self.regions_chargemult=regions_chargemult # List of list of charge,mult combos @@ -115,6 +127,11 @@ def __init__(self, theories_N=None, regions_N=None, regions_chargemult=None, for i,t in enumerate(self.theories_N): t.numcores=numcores print(f"Warning: Setting numcores={numcores} for Theory {i+1}: {t.theorynamelabel}") + #############################Dipanshu################## + if self.LL_full_theory is not None: + self.LL_full_theory.numcores=numcores + print(f"Warning: Setting numcores={numcores} for LL_full_theory: {self.LL_full_theory.theorynamelabel}") + ######################################################## else: print("Warning: numcores attribute was not set for ONIOMTheory.") print("This is fine, but check if numcores settings above are appropriate for each Theory object") @@ -419,9 +436,14 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non G_dict={} # (theory,region) -> gradient num_theories = len(self.theories_N) - # First doing LowLevel (LL) theory on Full region - ll_theory = self.theories_N[-1] - print(f"Running Theory LL ({ll_theory.theorynamelabel}) on Full-region ({len(full_elems)} atoms)") + ###########################Dipanshu############################ + # First doing LowLevel (LL) theory on Full region. + # Use LL_full_theory if the user provided one, otherwise fall back to theories_N[-1] (unchanged behavior) + if self.LL_full_theory is not None: + ll_theory = self.LL_full_theory + else: + ll_theory = self.theories_N[-1] + print(f"Running Theory LL_full ({ll_theory.theorynamelabel}) on Full-region ({len(full_elems)} atoms)") # Derive pointcharges unless full_pointcharges were already provided if self.full_pointcharges is None and self.embedding.lower() == "elstat": @@ -439,6 +461,8 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non ashexit() elif isinstance(ll_theory, xTBTheory): print(f"Theory is xTBTheory. Using default xtb charge model") + elif ll_theory.__class__.__name__ == "CP2KTheory": + print(f"Theory is CP2KTheory. Using {self.chargemodel} charge model") else: print("Problem: Theory-level not compatible with pointcharge-creation") ashexit() @@ -451,10 +475,14 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non # Changing filename here for GBW-file creation label = "LL_full" ll_theory_full.filename = f"{label}" - # RUN FULL + # RUN FULL (in its own subdirectory) + oniom_maindir = os.getcwd() + os.makedirs(label, exist_ok=True) + os.chdir(label) res_full = ll_theory_full.run(current_coords=full_coords, elems=full_elems, Grad=Grad, numcores=ll_theory.numcores, label=label, charge=self.fullregion_charge, mult=self.fullregion_mult) + os.chdir(oniom_maindir) if Grad: e_LL_full,g_LL_full = res_full else: @@ -474,6 +502,11 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non # Note: format issue # self.full_pointcharges = grabatomcharges_xTB_output(ll_theory.filename+'.out', chargemodel=self.chargemodel) self.full_pointcharges = grabatomcharges_xTB() + elif ll_theory_full.__class__.__name__ == "CP2KTheory": + from ash.interfaces.interface_CP2K import grabatomcharges_CP2K + cp2k_output = os.path.join( label, f"{ll_theory_full.filename}.out") + self.full_pointcharges = grabatomcharges_CP2K(cp2k_output, chargemodel=self.chargemodel) + print("self.full_pointcharges:", self.full_pointcharges) print(len(self.full_pointcharges)) @@ -630,6 +663,11 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non theory.filename = f"{label}" print(f"Running Theory {i+1} ({theory.theorynamelabel}) on Region {j+1} ({len(region_elems_final)} atoms)") + # Switch into a dedicated subdirectory for this theory/region combo + oniom_maindir = os.getcwd() + os.makedirs(label, exist_ok=True) + os.chdir(label) + # For an MM-theory like OpenMM we have to do some special handling if theory.theorytype == "MM": print("Case: Theory is MM") @@ -705,6 +743,8 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non res = theory.run(current_coords=region_coords_final, elems=region_elems_final, Grad=Grad, numcores=theory.numcores, PC=PC, current_MM_coords=pointchargecoords, MMcharges=pointcharges, mm_elems=mm_elems_for_qmprogram, label=label, charge=self.regions_chargemult[j][0], mult=self.regions_chargemult[j][1]) + os.chdir(oniom_maindir) + if PC and Grad: e,g,pg = res elif PC and not Grad: @@ -713,6 +753,7 @@ def run(self, current_coords=None, Grad=False, elems=None, charge=None, mult=Non e,g = res elif not PC and not Grad: e = res + # Keeping E and G in dicts E_dict[(i,j)] = e if Grad: