diff --git a/README.md b/README.md index 099789642..7b4eb651a 100644 --- a/README.md +++ b/README.md @@ -11,9 +11,8 @@ 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 `_ -R. Bjornsson*, J. Comput. Chem 2026, 47, e70359. +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. diff --git a/ash/__init__.py b/ash/__init__.py index ec90b8f75..a3ec98dd1 100644 --- a/ash/__init__.py +++ b/ash/__init__.py @@ -219,7 +219,7 @@ from ash.interfaces.interface_dlfind import DLFIND_optimizer,DLFIND_optimizerClass # Sella -from ash.interfaces.interface_sella import SellaOptimizer, SellaoptimizerClass +from ash.interfaces.interface_sella import SellaOptimizer, SellaoptimizerClass, SellaIRC # Other import ash.interfaces.interface_crest diff --git a/ash/dictionaries_lists.py b/ash/dictionaries_lists.py index 0996dc407..cea69daba 100644 --- a/ash/dictionaries_lists.py +++ b/ash/dictionaries_lists.py @@ -64,7 +64,7 @@ def __init__(self, name, symbol, atomnumber): self.atomnumber = atomnumber -element_dict_atname = {'h': Element('hydrogen', 'H', 1), 'he': Element('helium', 'He', 2), +element_dict_atname = {'m': Element('dummy', 'M', 0), 'h': Element('hydrogen', 'H', 1), 'he': Element('helium', 'He', 2), 'li': Element('lithium', 'Li', 3), 'be': Element('beryllium', 'Be', 4), 'b': Element('boron', 'B', 5), 'c': Element('carbon', 'C', 6), 'n': Element('nitrogen', 'N', 7), @@ -93,16 +93,31 @@ def __init__(self, name, symbol, atomnumber): 'sb': Element('antimony', 'Sb', 51), 'te': Element('tellurium', 'Te', 52), 'i': Element('iodine', 'I', 53), 'xe': Element('xenon', 'Xe', 54), 'cs': Element('cesium', 'Cs', 55), 'ba': Element('barium', 'Ba', 56), + 'la': Element('lanthanum', 'La', 57), 'ce': Element('cerium', 'Ce', 58), 'pr': Element('praseodymium', 'Pr', 59), + 'nd': Element('neodymium', 'Nd', 60), 'pm': Element('promethium', 'Pm', 61), 'sm': Element('samarium', 'Sm', 62), + 'eu': Element('europium', 'Eu', 63), 'gd': Element('gadolinium', 'Gd', 64), 'tb': Element('terbium', 'Tb', 65), + 'dy': Element('dysprosium', 'Dy', 66), 'ho': Element('holmium', 'Ho', 67), 'er': Element('erbium', 'Er', 68), + 'tm': Element('thulium', 'Tm', 69), 'yb': Element('ytterbium', 'Yb', 70), 'lu': Element('lutetium', 'Lu', 71), 'hf': Element('hafnium', 'Hf', 72), 'ta': Element('tantalum', 'Ta', 73), 'w': Element('tungsten', 'W', 74), 're': Element('rhenium', 'Re', 75), 'os': Element('osmium', 'Os', 76), 'ir': Element('iridium', 'Ir', 77), 'pt': Element('platinum', 'Pt', 78), 'au': Element('gold', 'Au', 79), 'hg': Element('mercury', 'Hg', 80), 'tl': Element('thallium', 'Tl', 81), 'pb': Element('lead', 'Pb', 82),'bi': Element('bismuth', 'Bi', 83),'po': Element('polonium', 'Po', 84), - 'at': Element('astatine', 'At', 85), 'rn': Element('radon', 'Rn', 86)} + 'at': Element('astatine', 'At', 85), 'rn': Element('radon', 'Rn', 86), 'fr': Element('francium', 'Fr', 87), + 'ra': Element('radium', 'Ra', 88), 'ac': Element('actinium', 'Ac', 89), + 'th': Element('thorium', 'Th', 90), 'pa': Element('protactinium', 'Pa', 91), 'u': Element('uranium', 'U', 92), + 'np': Element('neptunium', 'Np', 93), 'pu': Element('plutonium', 'Pu', 94), 'am': Element('americium', 'Am', 95), + 'cm': Element('curium', 'Cm', 96), 'bk': Element('berkelium', 'Bk', 97), 'cf': Element('californium', 'Cf', 98), + 'es': Element('einsteinium', 'Es', 99), 'fm': Element('fermium', 'Fm', 100), 'md': Element('mendelevium', 'Md', 101), + 'no': Element('nobelium', 'No', 102), 'lr': Element('lawrencium', 'Lr', 103), 'rf': Element('rutherfordium', 'Rf', 104), + 'db': Element('dubnium', 'Db', 105), 'sg': Element('seaborgium', 'Sg', 106), 'bh': Element('bohrium', 'Bh', 107), + 'hs': Element('hassium', 'Hs', 108), 'mt': Element('meitnerium', 'Mt', 109), 'ds': Element('darmstadtium', 'Ds', 110), + 'rg': Element('roentgenium', 'Rg', 111), 'cn': Element('copernicium', 'Cn', 112), 'nh': Element('nihonium', 'Nh', 113), + 'fl': Element('flerovium', 'Fl', 114), 'mc': Element('moscovium', 'Mc', 115), 'lv': Element('livermorium', 'Lv', 116), + 'ts': Element('tennessine', 'Ts', 117), 'og': Element('oganesson', 'Og', 118)} -# NOTE: redundant? remove? -element_dict_atnum = {1: Element('hydrogen', 'H', 1), 2: Element('helium', 'He', 2), 3: Element('lithium', 'Li', 3), +element_dict_atnum = {0: Element('dummy', 'M', 0), 1: Element('hydrogen', 'H', 1), 2: Element('helium', 'He', 2), 3: Element('lithium', 'Li', 3), 4: Element('beryllium', 'Be', 4), 5: Element('boron', 'B', 5), 6: Element('carbon', 'C', 6), 7: Element('nitrogen', 'N', 7), 8: Element('oxygen', 'O', 8), 9: Element('fluorine', 'F', 9), 10: Element('neon', 'Ne', 10), @@ -129,11 +144,31 @@ def __init__(self, name, symbol, atomnumber): 52: Element('tellurium', 'Te', 52), 53: Element('iodine', 'I', 53), 54: Element('xenon', 'Xe', 54), 55: Element('cesium', 'Cs', 55), 56: Element('barium', 'Ba', 56), + 57: Element('lanthanum', 'La', 57), 58: Element('cerium', 'Ce', 58), 59: Element('praseodymium', 'Pr', 59), + 60: Element('neodymium', 'Nd', 60), 61: Element('promethium', 'Pm', 61), 62: Element('samarium', 'Sm', 62), + 63: Element('europium', 'Eu', 63), 64: Element('gadolinium', 'Gd', 64), + 65: Element('terbium', 'Tb', 65), 66: Element('dysprosium', 'Dy', 66), + 67: Element('holmium', 'Ho', 67), 68: Element('erbium', 'Er', 68), + 69: Element('thulium', 'Tm', 69), 70: Element('ytterbium', 'Yb', 70), + 71: Element('lutetium', 'Lu', 71), 72: Element('hafnium', 'Hf', 72), 73: Element('tantalum', 'Ta', 73), 74: Element('tungsten', 'W', 74), 75: Element('rhenium', 'Re', 75), 76: Element('osmium', 'Os', 76), 77: Element('iridium', 'Ir', 77), 78: Element('platinum', 'Pt', 78), - 79: Element('gold', 'Au', 79), 80: Element('mercury', 'Hg', 80)} + 79: Element('gold', 'Au', 79), 80: Element('mercury', 'Hg', 80), + 81: Element('thallium', 'Tl', 81), 82: Element('lead', 'Pb', 82), + 83: Element('bismuth', 'Bi', 83), 84: Element('polonium', 'Po', 84), 85: Element('astatine', 'At', 85), + 86: Element('radon', 'Rn', 86), 87: Element('francium', 'Fr', 87), 88: Element('radium', 'Ra', 88), + 89: Element('actinium', 'Ac', 89), 90: Element('thorium', 'Th', 90), 91: Element('protactinium', 'Pa', 91), 92: Element('uranium', 'U', 92), + 93: Element('neptunium', 'Np', 93), 94: Element('plutonium', 'Pu', 94), 95: Element('americium', 'Am', 95), + 96: Element('curium', 'Cm', 96), 97: Element('berkelium', 'Bk', 97), 98: Element('californium', 'Cf', 98), + 99: Element('einsteinium', 'Es', 99), 100: Element('fermium', 'Fm', 100), 101: Element('mendelevium', 'Md', 101), + 102: Element('nobelium', 'No', 102), 103: Element('lawrencium', 'Lr', 103), 104: Element('rutherfordium', 'Rf', 104), + 105: Element('dubnium', 'Db', 105), 106: Element('seaborgium', 'Sg', 106), 107: Element('bohrium', 'Bh', 107), + 108: Element('hassium', 'Hs', 108), 109: Element('meitnerium', 'Mt', 109), 110: Element('darmstadtium', 'Ds', 110), + 111: Element('roentgenium', 'Rg', 111), 112: Element('copernicium', 'Cn', 112), 113: Element('nihonium', 'Nh', 113), + 114: Element('flerovium', 'Fl', 114), 115: Element('moscovium', 'Mc', 115), 116: Element('livermorium', 'Lv', 116), + 117: Element('tennessine', 'Ts', 117), 118: Element('oganesson', 'Og', 118)} # Ground-state spin multiplicities of all elements atom_spinmults = {'H': 2, 'He': 1, 'Li': 2, 'Be': 1, 'B': 2, 'C': 3, 'N': 4, 'O': 3, 'F': 2, 'Ne': 1, 'Na': 2, 'Mg': 1, diff --git a/ash/interfaces/interface_CFour.py b/ash/interfaces/interface_CFour.py index 5add5535a..908024704 100644 --- a/ash/interfaces/interface_CFour.py +++ b/ash/interfaces/interface_CFour.py @@ -151,7 +151,16 @@ def __init__(self, cfourdir=None, printlevel=2, cfouroptions=None, numcores=1, except shutil.SameFileError: pass - + #Storing absolute path of staged GENBAS file and caching ECPDATA contents in memory. + #Necessary because some job-types (e.g. NEB via Knarr, Singlepoint_parallel) change into + #subdirectories (image_N, worker dirs) before calling run(). The run-method will + #restore these files in the current working directory if they are missing there. + #Note: ECPDATA must be cached in memory since cleanup() (called below) deletes the file. + self.genbas_path = os.path.abspath('GENBAS') if os.path.isfile('GENBAS') else None + self.ecpdata_content = None + if os.path.isfile('ECPDATA'): + with open('ECPDATA') as ecpfile: + self.ecpdata_content = ecpfile.read() #Clean-up of possible old Cfour files before beginning #TODO: Skip cleanup of chosen files? @@ -321,6 +330,25 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el print("no current_coords") ashexit() + #Making sure GENBAS (and ECPDATA if available) is present in current working directory. + #Jobs like NEB (Knarr) and Singlepoint_parallel change into subdirectories (image_N, worker dirs) + #where the basis files staged by __init__ are not present, causing xcfour to crash. + if self.genbas_path is not None and not os.path.isfile("GENBAS"): + if os.path.isfile(self.genbas_path): + print(f"No GENBAS file found in current directory. Copying GENBAS from: {self.genbas_path}") + shutil.copyfile(self.genbas_path, "GENBAS") + if self.ecpdata_content is not None and not os.path.isfile("ECPDATA"): + with open("ECPDATA", "w") as ecpfile: + ecpfile.write(self.ecpdata_content) + + #Resetting EXTERN_POT for this run unless pointcharges are used (QM/MM job). + #Energy-only runs with PROP active set self.EXTERN_POT='ON' (dummy-pcharges trick to prevent + #reorientation), which is sticky on the object. A subsequent gradient run in a fresh directory + #(e.g. NEB image dirs) would then require a pcharges file that does not exist there, + #crashing CFour (mkvmol.f: End of file reading pcharges). + if PC is not True: + self.EXTERN_POT='OFF' + if self.DBOC is True: DBOC=True diff --git a/ash/interfaces/interface_ORCA.py b/ash/interfaces/interface_ORCA.py index bd7d2af06..9f9d9d387 100644 --- a/ash/interfaces/interface_ORCA.py +++ b/ash/interfaces/interface_ORCA.py @@ -1828,7 +1828,7 @@ def create_orca_input_pc(name,elems,coords,orcasimpleinput,orcablockinput,charge orcafile.write('! Freq' + '\n') if moreadfile is not None: print("MOREAD option active. Will read orbitals from file:", moreadfile) - orcafile.write('! MOREAD' + '\n') + orcafile.write('\n! MOREAD' + '\n') orcafile.write('%moinp \"{}\"'.format(moreadfile) + '\n') orcafile.write('%pointcharges "{}"\n'.format(pcfile)) orcafile.write(orcablockinput + '\n') @@ -1904,7 +1904,7 @@ def create_orca_input_plain(name,elems,coords,orcasimpleinput,orcablockinput,cha orcafile.write('! Freq' + '\n') if moreadfile is not None: print("MOREAD option active. Will read orbitals from file:", moreadfile) - orcafile.write('! MOREAD' + '\n') + orcafile.write('\n! MOREAD' + '\n') orcafile.write('%moinp \"{}\"'.format(moreadfile) + '\n') orcafile.write(orcablockinput) if atomstoflip is not None: @@ -2000,6 +2000,7 @@ def grabspinpop_ORCA(chargemodel,outputfile): coordgrab=False spinpops=[] BS=False #if broken-symmetry job + numatoms = int(pygrep('Number of atoms ...', outputfile)[-1]) #if if len(pygrep2("WARNING: Broken symmetry calculations", outputfile)): BS=True @@ -2032,7 +2033,9 @@ def grabspinpop_ORCA(chargemodel,outputfile): #If BS then we have grabbed charges for both high-spin and BS solution if BS is True: print("Broken-symmetry job detected. Only taking BS-state populations") - spinpops=spinpops[int(len(spinpops)/2):] + #spinpops=spinpops[int(len(spinpops)/2):] + if len(spinpops) != numatoms: + spinpops=spinpops[-numatoms:] #if len(spinpops) == 0: # print("Warning: No spinpopulations were found in ORCA outputfile") return spinpops diff --git a/ash/interfaces/interface_OpenMM.py b/ash/interfaces/interface_OpenMM.py index f62c99458..0ce9e112b 100644 --- a/ash/interfaces/interface_OpenMM.py +++ b/ash/interfaces/interface_OpenMM.py @@ -37,7 +37,7 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo Amberfiles=False, amberprmtopfile=None, properties=None, cluster_fragment=None, ASH_FF_file=None, nonbondedMethod_noPBC='NoCutoff', nonbonded_cutoff_noPBC=20, - xmlfiles=None, pdbfile=None, use_parmed=False, xmlsystemfile=None, + xmlfiles=None, pdbfile=None, pdbxfile=None, use_parmed=False, xmlsystemfile=None, do_energy_decomposition=False, periodic=False, periodic_cell_dimensions=None, PBCvectors=None, periodic_cell_vectors=None, charmm_periodic_cell_dimensions=None, customnonbondedforce=False, @@ -71,7 +71,6 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo self.analytic_hessian=False self.label=label self.fragment=fragment - # OPEN MM load try: import openmm @@ -241,7 +240,6 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo if periodic_cell_vectors is None: periodic_cell_vectors=PBCvectors - # #Always creates object we call self.forcefield that contains topology attribute if CHARMMfiles is True: if self.printlevel > 0: @@ -279,7 +277,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo self.atomtypes = [self.psf.atom_list[i].attype for i in range(0, len(self.psf.atom_list))] # TODO: Note: For atomnames it seems OpenMM converts atomnames to its own. Perhaps not useful self.atomnames = [self.psf.atom_list[i].name for i in range(0, len(self.psf.atom_list))] - self.mm_elements = [i.element.symbol for i in self.psf.topology.atoms()] + self.define_mm_elements(self.psf.topology) + #self.mm_elements = [i.element.symbol for i in self.psf.topology.atoms()] self.topology = self.psf.topology self.forcefield = self.psf @@ -322,8 +321,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo self.forcefield = self.grotop # TODO: Define resnames, resids, segmentnames, atomtypes, atomnames?? - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] - + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] elif Amberfiles is True: if self.printlevel > 0: print("Reading Amber files.") @@ -373,7 +372,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo #List of resids, resnames and mm_elements. Used by actregiondefine self.resids = [i.residue.index for i in self.prmtop.topology.atoms()] self.resnames = [i.residue.name for i in self.prmtop.topology.atoms()] - self.mm_elements = [i.element.symbol for i in self.prmtop.topology.atoms()] + self.define_mm_elements(self.prmtop.topology) + #self.mm_elements = [i.element.symbol for i in self.prmtop.topology.atoms()] self.atomnames = [i.name for i in self.prmtop.topology.atoms()] #NOTE: OpenMM does not grab Amber atomtypes for some reason. Feature request #TODO: Grab more topology information @@ -394,7 +394,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo # Check if PBC vectors in PDB-file pdb_pbc_vectors = pdb.topology.getPeriodicBoxVectors() self.forcefield = forcefield - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] elif ASH_FF_file is not None: if self.printlevel > 0: @@ -455,7 +456,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo residlabels=residlabels) pdb = openmm.app.PDBFile("cluster.pdb") self.topology = pdb.topology - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] self.forcefield = openmm.app.ForceField(xmlfile) # Load XMLfile for whole system @@ -483,7 +485,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo print("Reading topology from PDBfile:", pdbfile) pdb = openmm.app.PDBFile(pdbfile) self.topology = pdb.topology - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] # Check if PBC vectors in PDB-file pdb_pbc_vectors = pdb.topology.getPeriodicBoxVectors() # Simple OpenMM system without any forcefield defined. Requires ASH fragment @@ -509,21 +512,25 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo sigmas_per_res=[[0.0]*fragment.numatoms], epsilons_per_res=[[0.0]*fragment.numatoms], skip_nb=False) # Create dummy forcefield self.forcefield = openmm.app.ForceField(xmlfile) + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] - - # Read topology from PDB-file and XML-forcefield files to define forcefield + # Read topology from PDB-file or PDBx-file and XML-forcefield files to define forcefield else: if self.printlevel > 0: - print("Reading OpenMM XML forcefield files and PDB file") + print("Reading OpenMM XML forcefield files and PDB (or PDBx) file") print("xmlfiles:", str(xmlfiles).strip("[]")) print("pdbfile:", pdbfile) - if pdbfile is None: - print("Error:No pdbfile input provided") - ashexit() + print("pdbxfile:", pdbxfile) # This would be regular OpenMM Forcefield definition requiring XML file # Topology from PDBfile annoyingly enough - pdb = openmm.app.PDBFile(pdbfile) + if pdbfile is not None: + pdb = openmm.app.PDBFile(pdbfile) + elif pdbxfile is not None: + pdb = openmm.app.PDBxFile(pdbxfile) + else: + print("Error: No pdbfile or pdbxfile input provided") + ashexit() # Check if PBC vectors in PDB-file pdb_pbc_vectors = pdb.topology.getPeriodicBoxVectors() @@ -534,7 +541,8 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo self.resids = [i.residue.index for i in self.topology.atoms()] self.resnames = [i.residue.name for i in self.topology.atoms()] self.atomnames = [i.name for i in self.topology.atoms()] - self.mm_elements = [i.element.symbol for i in self.topology.atoms()] + self.define_mm_elements(self.topology) + #self.mm_elements = [i.element.symbol for i in self.topology.atoms()] # Dealing with possible user-defined residuetemplate_choice @@ -903,6 +911,15 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo print_time_rel(module_init_time, modulename="OpenMM object creation", moduleindex=3,currprintlevel=self.printlevel) + def define_mm_elements(self,topology): + try: + self.mm_elements = [i.element.symbol for i in topology.atoms()] + except: + print("Problem occurred while defining mm_elements.") + print("This may be due to virtual sites present") + print("mm_elements will be set to empty list") + self.mm_elements = [] + # Create a mixed MM/ML potential system from a ML potential (requires OpenMM-ML) # from openmmml import MLPotential # potential = MLPotential('ani2x') @@ -1635,6 +1652,14 @@ def create_integrator(self): self.integrator = openmm.VariableLangevinIntegrator(self.temperature * openmm.unit.kelvin, self.coupling_frequency / openmm.unit.picosecond, self.timestep * openmm.unit.picoseconds) + elif self.integrator_name == 'DrudeLangevinIntegrator': + print("here1") + # TODO: options + self.integrator = openmm.DrudeLangevinIntegrator(self.temperature * openmm.unit.kelvin, + self.coupling_frequency / openmm.unit.picosecond, + self.temperature * openmm.unit.kelvin, + self.timestep * openmm.unit.picoseconds,4) + print("here2") elif self.integrator_name == 'RPMDIntegrator': print("RPMDIntegrator will be used") print("Warning: Autoconstraints, rigidwater and other contraints must have been disabled.") @@ -2634,8 +2659,8 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi print("conda install -c conda-forge pdbfixer") ashexit() - if pdbfile == None: - print("You must provide a pdbfile= keyword argument") + if pdbfile is None: + print("You must provide a pdbfile keyword argument") ashexit() @@ -2981,6 +3006,7 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi fragment = Fragment(pdbfile="system_aftersolvent_ions.pdb") write_pdbfile_openMM(modeller.topology, modeller.positions, "finalsystem.pdb") + write_pdbxfile_openMM(modeller.topology, modeller.positions, "finalsystem.cif") fragment.print_system(filename="finalsystem.ygg") fragment.write_xyzfile(xyzfilename="finalsystem.xyz") @@ -3011,6 +3037,7 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi print("system_afterions.pdb and finalsystem.pdb (same)") print("\nFinal files:") print("finalsystem.pdb (PDB file)") + print("finalsystem.cif (PDBx/mmCIF file)") print("finalsystem.ygg (ASH fragment file)") print("finalsystem.xyz (XYZ coordinate file)") print("{} (System XML file)".format(systemxmlfile)) @@ -3018,15 +3045,20 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi print(BC.WARNING,"Strongly recommended: Check finalsystem.pdb carefully for correctness!", BC.END) print("\nTo use this system setup to define a future OpenMMTheory object you can either do:\n") - print(BC.OKMAGENTA,"1. Define using separate forcefield XML files:",BC.END) + print(BC.OKMAGENTA,"1. Define using separate forcefield XML files and PDB-file (for topology):",BC.END) if extraxmlfile is None: print(f"omm = OpenMMTheory(xmlfiles=[\"{xmlfile}\", \"{waterxmlfile}\"], pdbfile=\"finalsystem.pdb\", periodic={periodic})",BC.END) else: print(f"omm = OpenMMTheory(xmlfiles=[\"{xmlfile}\", \"{waterxmlfile}\", \"{extraxmlfile}\"], pdbfile=\"finalsystem.pdb\", periodic={periodic})",BC.END) - print(BC.OKMAGENTA,"2. Use forcefield object file :\n",BC.END, \ + print(BC.OKMAGENTA,"2. Define using separate forcefield XML files and PDBx/mmCIF file (instead of PDB):",BC.END) + if extraxmlfile is None: + print(f"omm = OpenMMTheory(xmlfiles=[\"{xmlfile}\", \"{waterxmlfile}\"], pdbxfile=\"finalsystem.cif\", periodic={periodic})",BC.END) + else: + print(f"omm = OpenMMTheory(xmlfiles=[\"{xmlfile}\", \"{waterxmlfile}\", \"{extraxmlfile}\"], pdbxfile=\"finalsystem.cif\", periodic={periodic})",BC.END) + print(BC.OKMAGENTA,"3. Use forcefield object file :\n",BC.END, \ f"omm = OpenMMTheory(topoforce=True, forcefield=forcefield_object, pdbfile=\"finalsystem.pdb\", topology=modeller.topology, periodic={periodic})",BC.END) - print(BC.OKMAGENTA,"3. Use full system XML-file (USUALLY NOT RECOMMENDED ):\n",BC.END, \ - f"omm = OpenMMTheory(xmlsystemfile=\"system_full.xml\", pdbfile=\"finalsystem.pdb\", periodic={periodic})",BC.END) + #print(BC.OKMAGENTA,"3. Use full system XML-file (USUALLY NOT RECOMMENDED ):\n",BC.END, \ + # f"omm = OpenMMTheory(xmlsystemfile=\"system_full.xml\", pdbfile=\"finalsystem.pdb\", periodic={periodic})",BC.END) print() print() if residuetemplate_choice is not None: @@ -3060,6 +3092,16 @@ def write_pdbfile_openMM(topology, positions, filename, connectivity_dict=None): openmm.app.PDBFile.writeFile(topology, positions, file=open(filename, 'w')) print("Wrote PDB-file:", filename) +def write_pdbxfile_openMM(topology, positions, filename, connectivity_dict=None): + import openmm.app + + if connectivity_dict is not None: + print("Connectivity passed to write_pdbxfile_openMM") + openmm_add_bonds_to_topology(topology,connectivity_dict) + + openmm.app.PDBxFile.writeFile(topology, positions, file=open(filename, 'w')) + print("Wrote PDBx-file:", filename) + #Take OpenMM topology and connectivity dictionary and add bonds to topology #in order for OpenMM PDBFile.writeFile to write CONECT lines def openmm_add_bonds_to_topology(topology,connectivity): @@ -4322,7 +4364,14 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m blastate.getPositions(asNumpy=True).value_in_unit( openmm.unit.angstrom), f) openmm.app.pdbfile.PDBFile.writeFooter(self.openmmobject.topology,f) - + # PDBx/mmCIF + pdbx_filename=self.trajfilename+"_firstframe.cif" + print("Writing intial frame to disk as PDBx/mmCIF-file:", pdbx_filename) + with open(pdbx_filename, 'w') as f: + openmm.app.pdbxfile.PDBxFile.writeHeader(self.openmmobject.topology, f) + openmm.app.pdbxfile.PDBxFile.writeModel(self.openmmobject.topology, + blastate.getPositions(asNumpy=True).value_in_unit( + openmm.unit.angstrom), f) ############################################################################### # MD LOOP for each Theory-Runtype: WRAP, QMMM, QM, ONIOM, dummy_MM, MM @@ -4879,9 +4928,9 @@ def finalize_simulation(self): #Topology (for header in PDB-files). Necessary self.openmmobject.topology.setPeriodicBoxVectors(self.state.getPeriodicBoxVectors()) - ######################################## - # Writing final frame to disk as PDB. - ######################################## + ################################################ + # Writing final frame to disk as PDB and PDBx + ################################################ pdb_filename=self.trajfilename+"_lastframe.pdb" print("Writing final frame to disk as PDB-file:", pdb_filename) with open(pdb_filename, 'w') as f: @@ -4891,7 +4940,17 @@ def finalize_simulation(self): openmm.unit.angstrom), f) openmm.app.pdbfile.PDBFile.writeFooter(self.openmmobject.topology,f) print(f"Trajectory : {self.trajfilename}.{self.trajectory_file_option}") - + # PDBx/mmCIF + pdbx_filename=self.trajfilename+"_lastframe.cif" + print("Writing final frame to disk as PDBx/mmCIF-file:", pdbx_filename) + with open(pdbx_filename, 'w') as f: + openmm.app.pdbxfile.PDBxFile.writeHeader(self.openmmobject.topology, f) + openmm.app.pdbxfile.PDBxFile.writeModel(self.openmmobject.topology, + self.state.getPositions(asNumpy=True).value_in_unit( + openmm.unit.angstrom), f) + print(f"Trajectory : {self.trajfilename}.{self.trajectory_file_option}") + + # Saving state to disk #Can be used to restart using statefile option print("Saving a statefile and checkpointfile of the final frame of the simulation: OpenMM_MD_final_state.xml and OpenMM_MD_final_checkpoint.chk") @@ -5045,8 +5104,8 @@ def OpenMM_box_equilibration(fragment=None, theory=None, datafilename="nptsim.cs MDtraj_imagetraj(f"{trajfilename}.dcd", f"{trajfilename}_lastframe.pdb") except ImportError: print("mdtraj library could not be imported. Skipping") - except ValueError: - print("mdtraj reimaging failed. Skipping") + except ValueError as e: + print(f"mdtraj reimaging failed. Skipping. Error: {e}") print_time_rel(module_init_time, modulename="OpenMM_box_equilibration", moduleindex=1) return md.state.getPeriodicBoxVectors() @@ -5510,8 +5569,8 @@ def Gentle_warm_up_MD(theory=None, fragment=None, time_steps=[0.0005,0.001,0.004 threshold=0.005, largest_values=10) except ImportError: print("mdtraj library could not be imported. Skipping") - except ValueError: - print("mdtraj reimaging failed. Skipping") + except ValueError as e: + print(f"mdtraj reimaging failed. Skipping. Error: {e}") print("Gentle_warm_up_MD finished successfully!") print_time_rel(module_init_time, modulename="Gentle_warm_up_MD", moduleindex=1) diff --git a/ash/interfaces/interface_crest.py b/ash/interfaces/interface_crest.py index 1cbee98a2..e6b396a1b 100644 --- a/ash/interfaces/interface_crest.py +++ b/ash/interfaces/interface_crest.py @@ -134,27 +134,166 @@ def new_call_crest(fragment=None, theory=None, crestdir=None, runtype="imtd-gc", pickle.dump(theory, open(theoryfilename, "wb" )) # Write ASH inputfile: ash_input.py - ashinput=f"""from ash import * -from ash.interfaces.interface_ORCA import print_gradient_in_ORCAformat + #ashinput=f"""from ash import * +#from ash.interfaces.interface_ORCA import print_gradient_in_ORCAformat +#import pickle +# +#frag = Fragment(xyzfile="genericinp.xyz", charge={charge},mult={mult}) +#Unpickling theory object +#theory = pickle.load(open(\"../{theoryfilename}\", \"rb\" )) +#result = Singlepoint(theory=theory, fragment=frag, Grad=True) +#print_gradient_in_ORCAformat(result.energy,result.gradient,"genericinp", extrabasename="") +#""" + ash_server_code = """# energy_server.py + +from ash import * import pickle +import socket +import json +import numpy as np -frag = Fragment(xyzfile="genericinp.xyz", charge={charge},mult={mult}) -#Unpickling theory object -theory = pickle.load(open(\"../{theoryfilename}\", \"rb\" )) -result = Singlepoint(theory=theory, fragment=frag, Grad=True) -print_gradient_in_ORCAformat(result.energy,result.gradient,"genericinp", extrabasename="") -""" - with open("ash_input.py", "w") as f: - f.write(ashinput) +print("Loading theory object...") - theorylines=f"""method = "generic" -binary = "python3 ../ash_input.py" -gradfile = "genericinp.engrad" -gradtype = "engrad" +theory = pickle.load(open("./theory.saved", "rb")) + +print("Theory loaded.") + +sock = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM) + +socket_path = "/tmp/ash_energy.sock" + +import os +if os.path.exists(socket_path): + os.remove(socket_path) + +sock.bind(socket_path) +sock.listen() + +print("Server ready") + +while True: + + conn, _ = sock.accept() + + try: + data = b"" + + while True: + chunk = conn.recv(4096) + + if not chunk: + break + + data += chunk + + request = json.loads(data.decode()) + + coords = np.array(request["coords"]) + + elements = request["elements"] + + frag = Fragment( + elems=elements, + coords=coords, + charge=0, + mult=1 + ) + + result = Singlepoint( + theory=theory, + fragment=frag, + Grad=True + ) + + response = { + "energy": float(result.energy), + "gradient": result.gradient.tolist() + } + + conn.sendall(json.dumps(response).encode()) + + except Exception as e: + + response = { + "error": str(e) + } + + conn.sendall(json.dumps(response).encode()) + + finally: + conn.close() + """ + + ash_client_code = """# energy_client.py +from ash import * +from ash.interfaces.interface_ORCA import print_gradient_in_ORCAformat + +import socket +import json + +frag = Fragment( + xyzfile="genericinp.xyz", + charge=0, + mult=1 +) + +request = { + "elements": frag.elems, + "coords": frag.coords.tolist() +} + +sock = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM) + +sock.connect("/tmp/ash_energy.sock") + +sock.sendall(json.dumps(request).encode()) + +sock.shutdown(socket.SHUT_WR) + +data = b"" + +while True: + + chunk = sock.recv(4096) + + if not chunk: + break + + data += chunk + +sock.close() + +response = json.loads(data.decode()) + +if "error" in response: + raise RuntimeError(response["error"]) + +energy = response["energy"] +gradient = response["gradient"] + +print_gradient_in_ORCAformat( + energy, + gradient, + "genericinp", + extrabasename="" +) """ + # Write server code + with open("ash_server.py", "w") as f: + f.write(ash_server_code) + + # Write client code + with open("ash_client.py", "w") as f: + f.write(ash_client_code) + # Write CREST toml file # Note: crest created dirs caleld calculation.level.X etc. and enters them + theorylines=f"""method = "generic" +binary = "python3 ../ash_client.py" +gradfile = "genericinp.engrad" +gradtype = "engrad" +""" tomlinput=f"""# CREST 3 input file input = "struc.xyz" runtype="{runtype}" @@ -185,15 +324,29 @@ def new_call_crest(fragment=None, theory=None, crestdir=None, runtype="imtd-gc", if runtype == "imtd-gc": print(f"Note:Energy window is {ewin} kcal/mol") - print("Now calling CREST like this: crest --input input.toml") - process = sp.run([crestdir + '/crest', '--input', 'input.toml']) + # Launching ASH server in background via subprocess + print("Launching ASH server in background...") + filex = open("ash_server.log", "w") + serverproc = sp.Popen(["python3", "ash_server.py"], stdout=filex, stderr=filex) + print("Now calling CREST like this: crest --input input.toml") + try: + process = sp.run([crestdir + '/crest', '--input', 'input.toml'], check=True) + finally: + print("Terminating ASH server...") + serverproc.terminate() + serverproc.wait(timeout=10) + filex.close() + print_time_rel(module_init_time, modulename='crest run', moduleindex=0) # Get conformers try: - list_conformers, list_energies = get_crest_conformers(charge=charge, mult=mult) - return list_conformers, list_energies + if runtype == "imtd-gc": + list_conformers, list_energies = get_crest_conformers(charge=charge, mult=mult) + return list_conformers, list_energies + else: + return None, None except: return None, None diff --git a/ash/interfaces/interface_dlfind.py b/ash/interfaces/interface_dlfind.py index a51f13342..febbfaa0a 100644 --- a/ash/interfaces/interface_dlfind.py +++ b/ash/interfaces/interface_dlfind.py @@ -101,34 +101,35 @@ def __init__(self,jobtype=None, fragment=None, fragment2=None, theory=None, char print("Error: You must either select a jobtype keyword (e.g. opt, neb, dimer, instanton) or select DL-FIND icoord and iopt codes") print("Example: DLFIND_optimizer(jobtype='opt') ") ashexit() - elif jobtype.lower() == "opt": - print("jobtype: opt chosen") - print("Choosing icoord=1 (HDLC internal coordinates) and iopt=3 (L-BFGS minimizer)") - print("For other coordinate-systems: choose icoord=0 (cartesian), icoord=2 (hdlc-tc), icoord=3 (dlc-prim), icoord=3 (dlc-tc)") - print("For other opt algorithms: choose iopt codes: 0: sd, 1: cg-autorestart, 2: cg-restart10, 3: lbfgs, 10: P-RFO") - icoord=1 - iopt=3 - elif jobtype.lower() == "tsopt" or jobtype.lower() == "ts": - print("jobtype: tsopt chosen") - print("Choosing icoord=3 (HDLC internal coordinates) and iopt=10 (P-RFO)") - print("Note: inithessian option is:", inithessian) - icoord=3 - iopt=10 - elif jobtype.lower() == "neb": - print("jobtype: neb chosen") - print("Choosing icoord=120 (NEB with frozen endpoints) and iopt=3 (L-BFGS)") - icoord=120 - iopt=3 - elif jobtype.lower() == "dimer": - print("jobtype: dimer chosen") - print("Choosing icoord=210 (Dimer) and iopt=3 (L-BFGS)") - icoord=210 - iopt=3 - elif jobtype.lower() == "qts" or jobtype.lower() == "instanton" : - print("jobtype: qts chosen (a.k.a. instanton)") - print("Choosing icoord=190 (qts) and iopt=3 (L-BFGS)") - icoord=190 - iopt=3 + elif jobtype is not None: + if jobtype.lower() == "opt": + print("jobtype: opt chosen") + print("Choosing icoord=1 (HDLC internal coordinates) and iopt=3 (L-BFGS minimizer)") + print("For other coordinate-systems: choose icoord=0 (cartesian), icoord=2 (hdlc-tc), icoord=3 (dlc-prim), icoord=3 (dlc-tc)") + print("For other opt algorithms: choose iopt codes: 0: sd, 1: cg-autorestart, 2: cg-restart10, 3: lbfgs, 10: P-RFO") + icoord=1 + iopt=3 + elif jobtype.lower() == "tsopt" or jobtype.lower() == "ts": + print("jobtype: tsopt chosen") + print("Choosing icoord=1 (HDLC internal coordinates) and iopt=10 (P-RFO)") + print("Note: inithessian option is:", inithessian) + icoord=1 + iopt=10 + elif jobtype.lower() == "neb": + print("jobtype: neb chosen") + print("Choosing icoord=120 (NEB with frozen endpoints) and iopt=3 (L-BFGS)") + icoord=120 + iopt=3 + elif jobtype.lower() == "dimer": + print("jobtype: dimer chosen") + print("Choosing icoord=210 (Dimer) and iopt=3 (L-BFGS)") + icoord=210 + iopt=3 + elif jobtype.lower() == "qts" or jobtype.lower() == "instanton" : + print("jobtype: qts chosen (a.k.a. instanton)") + print("Choosing icoord=190 (qts) and iopt=3 (L-BFGS)") + icoord=190 + iopt=3 else: print("No jobype selected.") print(f"Will start job based on chosen icoord={icoord} and iopt={iopt}") @@ -525,7 +526,7 @@ def setup_constraints_act_frozen(self): self.frozenatoms = self.frozenatoms+frozenatoms_xyz print("frozenatoms_z:", frozenatoms_z) if self.actatoms is not None: - print_if_level("Actatoms provided:", self.actatoms) + print_if_level("Actatoms provided: {self.actatoms}", self.printlevel,2) if self.PBC: print("PBC detected. Adding 4 dummy atoms to actatoms if not already present") diff --git a/ash/interfaces/interface_mace.py b/ash/interfaces/interface_mace.py index 56f55c8a1..ffa852533 100644 --- a/ash/interfaces/interface_mace.py +++ b/ash/interfaces/interface_mace.py @@ -429,14 +429,18 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el # Checking first f file exists self.check_file_exists(self.model_file) #Load model + print_time_rel(module_init_time, modulename=f'MACE run - before modelfile_load', moduleindex=2) self.modelfile_load() elif self.model_name is not None: print("Loading via model_name:", self.model_name) + print_time_rel(module_init_time, modulename=f'MACE run - before modelname_load', moduleindex=2) self.modelname_load() else: print("Error: Neither modelfile or modelname was defined.") ashexit() + print_time_rel(module_init_time, modulename=f'MACE run - after model-load', moduleindex=2) + # Creating ASE atoms object (MACE has ASE has dependency anyway) import ase if self.periodic: @@ -449,14 +453,17 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el # New simpler MACE interface via ASE # Works for foundational models + print_time_rel(module_init_time, modulename=f'MACE run - after atoms', moduleindex=2) if self.new_interface is True: # Add loaded model to ASE calculator atoms.calc = self.model # Run energy self.energy = atoms.get_potential_energy() * ash.constants.evtohar + print_time_rel(module_init_time, modulename=f'MACE run - after energy', moduleindex=2) print("Energy:", self.energy) forces = atoms.get_forces() + print_time_rel(module_init_time, modulename=f'MACE run - after forces', moduleindex=2) self.gradient = forces/-51.422067090480645 if self.stress: stress_ev_ang3 = atoms.get_stress(voigt=False) diff --git a/ash/interfaces/interface_pyscf.py b/ash/interfaces/interface_pyscf.py index c648e173f..10a4bb426 100644 --- a/ash/interfaces/interface_pyscf.py +++ b/ash/interfaces/interface_pyscf.py @@ -967,6 +967,12 @@ def read_chkfile(self,chkfile): chkfile_scftype="UHF" elif 2.0 in self.chkfileobject["mo_occ"]: chkfile_scftype="RHF" + else: + print("Occupations seems to fit neither RHF nor UHF") + print("Occupations:", self.chkfileobject["mo_occ"]) + print("Guessing GHF for now") + chkfile_scftype="GHF" + #Checking if mismatch between chkfile info and chosen scf-type #TODO: In principle we could convert RKS-info from chkfile to UKS-info and vice versa if chkfile_scftype == "UHF": @@ -977,7 +983,11 @@ def read_chkfile(self,chkfile): if self.scf_type == "UHF" or self.scf_type == "UKS": print("Warning: Mismatch between SCF-type in chkfile and PySCFTheory object. Ignoring chkfile") return False - + if chkfile_scftype == "GHF": + if self.scf_type not in ["GHF","GKS"]: + print("Warning: Mismatch between SCF-type in chkfile and PySCFTheory object. Ignoring chkfile") + return False + if chkfile_scftype == "UHF": #UNRESTRICTED if self.printlevel >= 1: @@ -1731,6 +1741,8 @@ def run_CC(self,mf=None, frozen_orbital_indices=None, CCmethod='CCSD(T)', CC_dir print("Now starting CCSD calculation") if self.scf_type == "RHF": cc = pyscf_cc.CCSD(mf, frozen_orbital_indices,mo_coeff=mo_coefficients) + elif self.scf_type == "GHF": + cc = pyscf_cc.GCCSD(mf, frozen_orbital_indices,mo_coeff=mo_coefficients) elif self.scf_type == "ROHF": cc = pyscf_cc.CCSD(mf, frozen_orbital_indices,mo_coeff=mo_coefficients) elif self.scf_type == "UHF": @@ -2082,12 +2094,7 @@ def create_mf(self): #RKS v UKS v RHF v UHF v GHF v GKS #TODO: Dirac HF and KS also - if "skala" in self.functional.lower(): - print("here") - from skala.pyscf import SkalaKS - self.mf = SkalaKS(self.molcellobject, xc=self.functional) - - elif self.scf_type == 'RKS': + if self.scf_type == 'RKS': self.mf = scf.RKS(self.molcellobject) elif self.scf_type == 'ROKS': self.mf = scf.ROKS(self.molcellobject) @@ -2137,6 +2144,13 @@ def create_mf(self): else: print("Unknown scf-type:", self.scf_type) ashexit() + + # If Skala then override + if self.functional is not None and "skala" in self.functional.lower(): + from skala.pyscf import SkalaKS + self.mf = SkalaKS(self.molcellobject, xc=self.functional) + + print("mf object:", self.mf) #Probably depreceated. Created mf for GPU. @@ -2528,6 +2542,7 @@ def run_BS_SCF(self, mult=None, dm=None): #Independent method to run SCF using previously defined mf object and possible input dm def run_SCF(self,mf=None, dm=None, max_cycle=None): import pyscf + import pyscf.dft if self.printlevel >= 1: print("\nInside run_SCF") module_init_time=time.time() @@ -2588,9 +2603,10 @@ def run_SCF(self,mf=None, dm=None, max_cycle=None): import pyscf.pbc if isinstance(self.mf, pyscf.pbc.dft.rks.RKS): self.num_orbs = len(self.mf.mo_occ) # Restricted + elif self.scf_type in ['UKS','UHF']: + self.num_orbs = len(self.mf.mo_occ[0]) # Unrestricted else: - #UHF/UKS - self.num_orbs = len(self.mf.mo_occ[0]) + self.num_orbs = len(self.mf.mo_occ) # GHF/GKS if self.printlevel >= 1: print("Number of orbitals:", self.num_orbs) @@ -2687,7 +2703,6 @@ def prepare_run(self, current_coords=None, current_MM_coords=None, MMcharges=Non if self.printlevel >= 1: print("Number of electrons:", self.num_electrons) print() - ############################### #CREATE MOL OBJECT or CELL ############################### @@ -2767,7 +2782,13 @@ def prepare_run(self, current_coords=None, current_MM_coords=None, MMcharges=Non #Convert non-relativistic mf object to spin-free x2c if self.x2c is True if self.x2c is True: print("x2c is True. Changing SCF object to relativistic x2c Hamiltonian") - self.mf = self.mf.sfx2c1e() + if self.scf_type in ['RHF','UHF','RKS','UKS']: + print("Warning:Scalar relativistic x2c (spin-free) will be used") + self.mf = self.mf.sfx2c1e() + elif self.scf_type in ['GHF','GKS']: + print("Warning: x2c with spin orbit coupling enabled") + self.mf = self.mf.x2c() + ########### # PRINTING @@ -2918,7 +2939,8 @@ def actualrun(self, current_coords=None, current_MM_coords=None, MMcharges=None, if self.printlevel >1: print("Total num. orbitals:", self.num_scf_orbitals_alpha) if self.printlevel >1: - self.run_population_analysis(self.mf, dm=None, unrestricted=False, type='Mulliken', label='SCF') + if self.do_pop_analysis: + self.run_population_analysis(self.mf, dm=None, unrestricted=False, type='Mulliken', label='SCF') else: #UHF/UKS self.num_scf_orbitals_alpha=len(scf_result.mo_occ[0]) @@ -3050,9 +3072,22 @@ def actualrun(self, current_coords=None, current_MM_coords=None, MMcharges=None, if self.printlevel >1: print("MOM-SCF Gradient calculation done") else: - print("Gradient for postSCF methods is not implemented in ASH interface") - #TODO: Enable TDDFT, CASSCF, MP2, CC gradient etc - ashexit() + print("Warning: Gradient for postSCF methods in the pySCF interface are currently experimental.") + if self.CC: + if self.CCmethod == "CCSD(T)": + from pyscf.grad import ccsd_t as ccsd_t_grad + g = ccsd_t_grad.Gradients(self.ccobject).kernel() + print('CCSD(T) nuclear gradient:', g) + elif self.CCmethod == "CCSD": + from pyscf.grad import ccsd as ccsd_grad + g = ccsd_grad.Gradients(self.ccobject).kernel() + print('CCSD nuclear gradient:', g) + else: + print("CC method not recognized for gradient calculation. No gradient calculated.") + else: + print("No post-SCF method recognized for gradient calculation. No gradient calculated.") + self.gradient=None + #Calculate regular SCF gradient else: if self.printlevel >1: @@ -3475,8 +3510,8 @@ def KS_inversion_n2v(pyscftheoryobj, dm, method='PDECO', numcores=1, opt_max_ite import n2v import gbasis import pylibxc2 - except ModuleNotFoundError: - print("ModuleNotFoundError:") + except ModuleNotFoundError as e: + print("ModuleNotFoundError:", e) print("KS_inversion_n2v requires installation of n2v module and additional packages") print("See https://github.com/wasserman-group/n2v for details") print("""\n#Install pylibxc2 @@ -3533,7 +3568,7 @@ def KS_inversion_n2v(pyscftheoryobj, dm, method='PDECO', numcores=1, opt_max_ite #Takes pyscfheoryobject and DM as input, solves the inversion problem and returns MO coefficients, occupations,energies and new DM def KS_inversion_kspies(pyscftheoryobj, dm, numcores=1, method='WY', WY_method='trust-exact', pbas=None, - ZMP_lambda=128, ZMP_levelshift=True, ZMP_LS_scaling=8, + ZMP_lambda=128, ZMP_levelshift=True, ZMP_LS_scaling=0.1, ZMP_cycles=400, guide='faxc', DF=True, vxc_method = 'eval_vh', plot_vxc=False, vxc_coords=None, chosen_axes=None, @@ -3670,8 +3705,8 @@ def run_eval_vh(pyscftheoryobj, kspiesobj,coords_bohr, dm,l=0): for l in [ ZMP_lambda/8, ZMP_lambda/4, ZMP_lambda/2, ZMP_lambda]: #for l in [ZMP_lambda/(ZMP_LS_scaling*ZMP_LS_scaling), ZMP_lambda/ZMP_LS_scaling,ZMP_lambda]: print("Lambda:", l) - print("Levelshift (lambda*0.1):", l*0.1) - zmp_a.level_shift = l*0.1 + print("Levelshift (lambda*levelshiftscaling):", l*ZMP_LS_scaling) + zmp_a.level_shift = l*ZMP_LS_scaling zmp_a.zscf(l) if plot_vxc is True and plot_all_lambdas is True: diff --git a/ash/interfaces/interface_sella.py b/ash/interfaces/interface_sella.py index c0ff86b0d..150a5f5a5 100644 --- a/ash/interfaces/interface_sella.py +++ b/ash/interfaces/interface_sella.py @@ -12,8 +12,33 @@ from ash.modules.module_QMMM import QMMMTheory # Sella TS optimizer -# TODO active region -# TODO PBC + +def SellaIRC(theory=None, fragment=None, charge=None, mult=None, printlevel=2, NumGrad=False, + convergence_gmax=1e-4, maxiter=150, result_write_to_disk=False, + constraints=None, actatoms=None, frozenatoms=None, + gamma=0.03, eta=1e-4, IRC_dx=0.1): + """ + Wrapper function around SellaoptimizerClass + """ + timeA=time.time() + + optimizer = SellaoptimizerClass(convergence_gmax=convergence_gmax, + printlevel=printlevel, maxiter=maxiter, result_write_to_disk=result_write_to_disk, + constraints=constraints, actatoms=actatoms, frozenatoms=frozenatoms, + gamma=gamma, eta=eta, IRC_dx=IRC_dx) + + # If NumGrad then we wrap theory object into NumGrad class object + if NumGrad: + print("NumGrad flag detected. Wrapping theory object into NumGrad class") + print("This enables numerical-gradient calculation for theory") + theory = NumGradclass(theory=theory) + + # Providing theory and fragment to run method. + result = optimizer.IRCrun(theory=theory, fragment=fragment, charge=charge, mult=mult) + if printlevel >= 1: + print_time_rel(timeA, modulename='Sella', moduleindex=1) + + return result def SellaOptimizer(theory=None, fragment=None, charge=None, mult=None, printlevel=2, NumGrad=False, @@ -53,7 +78,7 @@ class SellaoptimizerClass: def __init__(self,printlevel=2, convergence_gmax=3e-4, maxiter=150, result_write_to_disk=False, constraints=None, actatoms=None, frozenatoms=None, - gamma=0.03, eta=1e-4): + gamma=0.03, eta=1e-4, IRC_dx=0.1): self.printlevel=printlevel print_line_with_mainheader("SellaOptimizer initialization") @@ -77,7 +102,7 @@ def __init__(self,printlevel=2, self.actatoms = actatoms self.gamma = gamma self.eta = eta - + self.IRC_dx = IRC_dx print_if_level(f"GradMax convergence tolerance: {self.convergence_gmax} Eh/Bohr", self.printlevel, 2) print_if_level(f"Converted tolerance for Sella: {self.tolerance_ev_ang} eV/Angstrom", self.printlevel, 2) @@ -159,6 +184,7 @@ def setup_constraints(self, atoms, constraints, fragment): print("All Sella constraints:", sellacons) return sellacons + def run(self, theory=None, fragment=None, charge=None, mult=None,constraints=None): print_line_with_subheader1("Running Sella optimization") @@ -201,6 +227,7 @@ def run(self, theory=None, fragment=None, charge=None, mult=None,constraints=Non atoms, constraints=sella_constraints, gamma=self.gamma, eta=self.eta) + def write_traj(a=atoms, trajname="sella_optim"): print(f"Writing (active) trajectory to file: {trajname}.xyz") self.active_fragment.coords = copy.copy(a.get_positions()) @@ -217,7 +244,6 @@ def write_qmregion_traj(a=atoms, trajname="sella_optim_qmregion"): frag = Fragment(coords=qm_coords, elems=qm_elems, printlevel=0) frag.write_xyzfile(xyzfilename=trajname+'.xyz', writemode='a') - # Attaching traj function #dyn.attach(print_step, interval=1) dyn.attach(write_traj, interval=1) @@ -286,8 +312,101 @@ def write_qmregion_traj(a=atoms, trajname="sella_optim_qmregion"): if self.result_write_to_disk is True: result.write_to_disk(filename="SellaOptimizer.result") return result + + def IRCrun(self, theory=None, fragment=None, charge=None, mult=None,constraints=None): + + print_line_with_subheader1("Running Sella IRC") + from sella import IRC + import ase + + # Active region setup. For a big system, we have to pass only the active region geometry to Sella + if self.actatoms is not None: + self.original_fragment = copy.deepcopy(fragment) + self.active_fragment = self.setup_active_region_geometry(fragment) + print(f"Active region fragment contains {self.active_fragment.numatoms} atoms") + else: + self.original_fragment=None # + self.active_fragment = fragment + # Creating ASE object + fragment.printlevel=0 + atoms = ase.atoms.Atoms(self.active_fragment.elems,positions=self.active_fragment.coords) + + # Setup constraints for Sella + #sella_constraints = None + #if self.constraints is not None or self.frozenatoms is not None: + # sella_constraints = self.setup_constraints(atoms, constraints,fragment) + #print("sella_constraints:", sella_constraints) + + # Attaching calculator + print("Creating ASH-ASE calculator") + atoms.calc = ASH_ASE_calculator(theory=theory, fragment=self.active_fragment, + full_fragment=self.original_fragment, actatoms=self.actatoms) + + # Attaching traj functions + def write_traj(a=atoms, trajname="sella_IRC"): + print(f"Writing (active) trajectory to file: {trajname}.xyz") + self.active_fragment.coords = copy.copy(a.get_positions()) + self.active_fragment.write_xyzfile(xyzfilename=trajname+'.xyz', writemode='a') + def write_full_traj(a=atoms, trajname="sella_optim_full"): + print(f"Writing full trajectory to file: {trajname}.xyz") + #self.original_fragment = copy.copy(a.get_positions()) + atoms.calc.full_fragment.write_xyzfile(xyzfilename=trajname+'.xyz', writemode='a') + def write_qmregion_traj(a=atoms, trajname="sella_optim_qmregion"): + print(f"Writing QM-region trajectory to file: {trajname}.xyz") + qm_elems = [atoms.calc.full_fragment.elems[i] for i in theory.qmatoms] + qm_coords = np.array([atoms.calc.full_fragment.coords[i] for i in theory.qmatoms]) + frag = Fragment(coords=qm_coords, elems=qm_elems, printlevel=0) + frag.write_xyzfile(xyzfilename=trajname+'.xyz', writemode='a') + + print("Creating Sella IRC object.") + print("dx:", self.IRC_dx, "eta:", self.eta, "gamma:", self.gamma) + opt = IRC(atoms, trajectory='irc.traj', dx=self.IRC_dx, ninner_iter=10, + eta=self.eta, gamma=self.gamma) + # Attaching traj function + opt.attach(write_traj, interval=1) + # Attaching full traj write also if using active region + if self.actatoms is not None: + opt.attach(write_full_traj, interval=1) + if isinstance(theory, QMMMTheory): + opt.attach(write_qmregion_traj, interval=1) + + print("Now running Sella IRC in forward direction.") + print(f"Convergence fmax: {self.tolerance_ev_ang} eV/Å, steps: {self.maxiter}") + # Running forward IRC step by step + conv = False + try: + conv = opt.run(fmax=self.tolerance_ev_ang, steps=self.maxiter, direction='forward') + except Exception as e: + print("Error occurred while running forward IRC:",e) + print("Continuing...") + print("Forward IRC run completed") + if conv: + print("Forward IRC converged successfully") + else: + print("Warning: Forward IRC did not converge. Continuing nonetheless.") + + print("Now running Sella IRC in reverse direction.") + print(f"Convergence fmax: {self.tolerance_ev_ang} eV/Å, steps: {self.maxiter}") + conv=False + try: + conv = opt.run(fmax=self.tolerance_ev_ang, steps=self.maxiter, direction='reverse') + except Exception as e: + print("Error occurred while running reverse IRC:",e) + print("Continuing...") + print("Reverse IRC run completed.") + if conv: + print("Reverse IRC converged successfully") + else: + print("Warning: Reverse IRC did not converge.") + + + result = ASH_Results(label="SellaIRC") + if self.result_write_to_disk is True: + result.write_to_disk(filename="SellaIRC.result") + return result + # Simpler ASH-ASE calculator class ASH_ASE_calculator: def __init__(self, theory=None, fragment=None, full_fragment=None, actatoms=None): @@ -305,9 +424,9 @@ def __init__(self, theory=None, fragment=None, full_fragment=None, actatoms=None # Initializing coordinates used by Sella self.coords=fragment.coords - def get_potential_energy(self, atomsobj): - #print("Called ASHcalc get_potential_energy") - #print("Energy call number:", self.energycalls) + def get_potential_energy(self, atomsobj, force_consistent=None): + print("Called ASHcalc get_potential_energy") + print("Energy call number:", self.energycalls) self.energycalls += 1 # Have coordinates changed? if np.array_equal(atomsobj.get_positions(), self.coords): @@ -316,8 +435,7 @@ def get_potential_energy(self, atomsobj): return self.energy_eV else: print("No energy available (1st step?). Will do calculation") - # ? - exit() + self.get_forces(atomsobj) return self.energy_eV def get_forces(self, atomsobj): @@ -361,5 +479,7 @@ def get_forces(self, atomsobj): # Energy self.energy_eH = energy self.energy_eV = energy*hartoeV - + #print("Energy:", self.energy_eH) + #print("Forces:", self.forces) return self.forces + diff --git a/ash/interfaces/interface_xtb.py b/ash/interfaces/interface_xtb.py index 84aeb18c3..05151cd5f 100644 --- a/ash/interfaces/interface_xtb.py +++ b/ash/interfaces/interface_xtb.py @@ -452,7 +452,7 @@ def get_cell_gradient(self): def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None, mm_elems=None, printlevel=None, elems=None, Grad=False, PC=False, numcores=None, label=None, charge=None, mult=None): module_init_time=time.time() - + if self.runmode == 'library': from xtb.interface import Calculator, Param diff --git a/ash/modules/module_PES_rewrite.py b/ash/modules/module_PES_rewrite.py index 447327b9d..8caa1470f 100644 --- a/ash/modules/module_PES_rewrite.py +++ b/ash/modules/module_PES_rewrite.py @@ -41,7 +41,7 @@ def PhotoElectron(theory=None, fragment=None, method=None, vibrational_option=No MRCI_CASCI_Final=False, MRCI_SOC=False, btPNO=False, DLPNO=False, no_shakeup=False, virt_offset=0, path_wfoverlap=None, tprintwfvalue=1e-5, noDyson=False, - OODFT_CC=False): + OODFT_CC=False, deltascf_moread_init=True): """ Wrapper function around PhotoElectron Class """ @@ -56,10 +56,10 @@ def PhotoElectron(theory=None, fragment=None, method=None, vibrational_option=No Ionizedstate_charge=Ionizedstate_charge, Ionizedstate_mult=Ionizedstate_mult, numionstates=numionstates, initialorbitalfiles=initialorbitalfiles, densities=densities, densgridvalue=densgridvalue, tda=tda,brokensym=brokensym, HSmult=HSmult, atomstoflip=atomstoflip, check_stability=check_stability, - deltaSCF_ionize=deltaSCF_ionize, deltaSCF_PMOM=deltaSCF_ionize, deltaSCFkeyword=deltaSCFkeyword, + deltaSCF_ionize=deltaSCF_ionize, deltaSCF_PMOM=deltaSCF_PMOM, deltaSCFkeyword=deltaSCFkeyword, CAS_Initial=CAS_Initial, CAS_Final=CAS_Final, no_shakeup=no_shakeup,virt_offset=virt_offset, MRCI_CASCI_Final=MRCI_CASCI_Final, MRCI_SOC=MRCI_SOC, CASCI_Final=CASCI_Final, - btPNO=btPNO, DLPNO=DLPNO, + btPNO=btPNO, DLPNO=DLPNO, deltascf_moread_init=deltascf_moread_init, path_wfoverlap=path_wfoverlap, tprintwfvalue=tprintwfvalue, noDyson=noDyson, OODFT_CC=OODFT_CC) result = photo.run() @@ -79,7 +79,7 @@ def __init__(self,theory=None, fragment=None, method=None, vibrational_option=No MRCI_CASCI_Final=False, MRCI_SOC=False, CASCI_Final=False, btPNO=False, DLPNO=False, path_wfoverlap=None, tprintwfvalue=1e-5, noDyson=False, - OODFT_CC=False): + OODFT_CC=False, deltascf_moread_init=True): """ PhotoElectron module """ @@ -201,6 +201,8 @@ def __init__(self,theory=None, fragment=None, method=None, vibrational_option=No self.deltaSCF_PMOM=deltaSCF_PMOM #Whether to use PMOM or not self.deltaSCFkeyword=deltaSCFkeyword # Add extra ORCA simple keyword when doing deltaSCF calcs only self.OODFT_CC=OODFT_CC # CCSD(T) on top of deltaSCF + self.deltascf_moread_init = deltascf_moread_init # Whether to add MOREAD line to read initial state GBW file for each ionized state in deltaSCF calcs + print("PES method:", self.method) if self.method == 'MRCI' or self.method=='MREOM': print("MREOM:", self.MREOM) @@ -348,15 +350,18 @@ def setup_ORCA_object(self): self.theory.orcasimpleinput = self.theory.orcasimpleinput + ' Normalprint' if self.method == "OODFT": - if 'UKS' not in self.theory.orcasimpleinput.upper(): - self.theory.orcasimpleinput = self.theory.orcasimpleinput + ' UKS' + if 'UHF' not in self.theory.orcasimpleinput.upper(): + # Avoid adding UHF if UHF or HF present + #if 'UHF' not in self.theory.orcasimpleinput.upper() and 'HF' not in self.theory.orcasimpleinput.upper(): + self.theory.orcasimpleinput = self.theory.orcasimpleinput + ' UHF' + if self.brokensym is True: self.theory.brokensym=True self.theory.HSmult=self.HSmult self.theory.atomstoflip=self.atomstoflip - #Making sure UKS always present if brokensym feature active. Important for open-shell singlets - if 'UKS' not in self.theory.orcasimpleinput.upper(): - self.theory.orcasimpleinput = self.theory.orcasimpleinput + ' UKS' + #Making sure UHF always present if brokensym feature active. Important for open-shell singlets + if 'UHF' not in self.theory.orcasimpleinput.upper(): + self.theory.orcasimpleinput = self.theory.orcasimpleinput + ' UHF' #Preserve original self.orig_orcasimpleinput=copy.copy(self.theory.orcasimpleinput) @@ -1111,6 +1116,7 @@ def run_DELTASCF(self,fragment,theory): print("now done with initial state SCF") + #Create strings for the SCF configurations and deltaSCF lines if self.deltaSCF_ionize is True: print("deltaSCF_ionize option active!") @@ -1124,13 +1130,12 @@ def run_DELTASCF(self,fragment,theory): IPs_all=[] Ionstates_energies_all=[] - # DELTASCF extra keywords - if self.deltaSCFkeyword is not None: - print("Adding deltaSCFkeyword to ORCA input string:", self.deltaSCFkeyword) - theory.orcasimpleinput += f" {self.deltaSCFkeyword} " + print("self.Finalstates:", self.Finalstates) #LOOPING over Finalstate-multiplicities for fstate in self.Finalstates: + print("fstate.mult:", fstate.mult) + if self.deltaSCF_ionize: print("deltaSCF_ionize option active. Note: this uses initial-state charge/multiplicity in ORCA input") charge = self.Initialstate_charge @@ -1143,28 +1148,72 @@ def run_DELTASCF(self,fragment,theory): print(f"\nNow running Finalstate MULT {fstate.mult} state-number {i}") label=f"{fstate.label}_state{i}" print("Label:", label) + + # DELTASCF extra keywords + if self.deltaSCFkeyword is not None: + for kw in self.deltaSCFkeyword.split(): + if kw not in theory.orcasimpleinput: + print(f"Adding {kw} to ORCA input string for DELTASCF calculation.") + # Skipping GMF if first ionized state (ORCA bug) + if kw.upper() == "GMF" and i == 0: + print("deltaSCFkeyword included GMF and this is the first ionized state. Skipping GMF for this state due to ORCA bug.") + else: + theory.orcasimpleinput += f" {kw} " + else: + print(f"{kw} already in ORCA input string. Not adding again.") + #Creating DELTASCF block if fstate.mult > self.stateI.mult: deltascfblock=f"%SCF \n PMOM {self.deltaSCF_PMOM} \n{deltascfline_CFG_betahole[i]}\nEND" else: deltascfblock=f"%SCF \n PMOM {self.deltaSCF_PMOM} \n{deltascfline_CFG_alphahole[i]}\nEND" - # Adding UKS if not - if 'UKS' not in theory.orcasimpleinput: - theory.orcasimpleinput += " UKS " + # Adding UHF if not + if 'UHF' not in theory.orcasimpleinput and 'HF' not in theory.orcasimpleinput: + theory.orcasimpleinput += " UHF " #Adding DELTASCF keyword to simpleinput (fine even for ground-state) if 'DELTASCF' not in theory.extraline: theory.extraline="! DELTASCF" #Adding ALPHACONF/BETACONF line as separate SCF block (empty if ground ion state) theory.orcablocks = self.orig_orcablocks + deltascfblock + + # MOREAD INIT state + if self.deltascf_moread_init is True: + print("deltascf_moread_init is True. Adding MOREAD line to read initial state GBW file for each ionized state.") + theory.moreadfile="Init_State.gbw" + state_result = ash.Singlepoint(fragment=fragment, theory=theory, charge=charge, mult=mult) + + #Keeping copy of input/outputfile and GBW file + shutil.copyfile(theory.filename + '.out', './' + label + '.out') + shutil.copyfile(theory.filename + '.inp', './' + label + '.inp') + shutil.copyfile(theory.filename + '.gbw', './' + label + '.gbw') + # CCSD(T) correction on top print("self.OODFT_CC:", self.OODFT_CC) if self.OODFT_CC: print("Now running noiter CCSD(T) on top of deltaSCF") theory.extraline = theory.extraline.replace("DELTASCF","CCSD(T) noiter ") + theory.orcasimpleinput = theory.orcasimpleinput.replace("FRSOSCF","") + theory.orcasimpleinput = theory.orcasimpleinput.replace("FreezeAndRelease","") + + + # Removing SCF maxiter line if present in theory.orcablocks to avoid conflict with CCSD(T) maxiter line + orig_orcablocks = self.orig_orcablocks + if "maxiter" in theory.orcablocks: + print("Removing maxiter line from ORCA blocks to avoid conflict with CCSD(T) maxiter") + orig_orcablocks = '\n'.join([line for line in theory.orcablocks.splitlines() if "maxiter" not in line]) + + # mdci block + theory.orcablocks = orig_orcablocks + f"""\n%mdci +maxiter 100 +maxdiis 20 +end +""" + state_result = ash.Singlepoint(fragment=fragment, theory=theory, charge=charge, mult=mult) + shutil.copyfile(theory.filename + '.out', './' + label + 'cc_noiter.out') finalsinglepointenergy = state_result.energy ip = (finalsinglepointenergy-self.stateI.energy)*ash.constants.hartoeV @@ -1184,10 +1233,7 @@ def run_DELTASCF(self,fragment,theory): #print("Initial state occupied MO energies (alpha):", self.stateI.occorbs_alpha) #print("Initial state SCF-type:", self.stateI.hftyp) - #Keeping copy of input/outputfile and GBW file - shutil.copyfile(theory.filename + '.out', './' + label + '.out') - shutil.copyfile(theory.filename + '.inp', './' + label + '.inp') - shutil.copyfile(theory.filename + '.gbw', './' + label + '.gbw') + #Calculate SCF eldensity and spindensity if requested if self.densities == 'SCF' or self.densities == 'All': @@ -1664,10 +1710,13 @@ def run_SCF_InitState(self,fragment,theory): if self.OODFT_CC: print("SCF InitState. Now running noiter CCSD(T) on top of deltaSCF") + # Removing SCF maxiter line if present in theory.orcablocks to avoid conflict with CCSD(T) maxiter line + if "maxiter" in theory.orcablocks: + print("Removing maxiter line from ORCA blocks to avoid conflict with CCSD(T) maxiter") + theory.orcablocks = '\n'.join([line for line in theory.orcablocks.splitlines() if "maxiter" not in line]) theory.extraline = theory.extraline + "! CCSD(T) noiter " state_result = ash.Singlepoint(fragment=fragment, theory=theory, charge=self.Initialstate_charge, mult=self.Initialstate_mult) - #Grab energy of initial state if self.method == 'CASSCF' or self.method =='CASCI': self.stateI.energy=casscfenergygrab(theory.filename+'.out') @@ -1680,6 +1729,7 @@ def run_SCF_InitState(self,fragment,theory): self.stateI.occorbs_alpha, self.stateI.occorbs_beta, self.stateI.hftyp = orbitalgrab(theory.filename+'.out') print("Initial state occupied MO energies (alpha):", self.stateI.occorbs_alpha) + print("Initial state occupied MO energies (beta):", self.stateI.occorbs_beta) print("Initial state SCF-type:", self.stateI.hftyp) # Specify whether Initial state is restricted or not. if self.stateI.hftyp == "UHF": diff --git a/ash/modules/module_QMMM.py b/ash/modules/module_QMMM.py index b85631b72..04cf057a1 100644 --- a/ash/modules/module_QMMM.py +++ b/ash/modules/module_QMMM.py @@ -176,6 +176,9 @@ def __init__(self, qm_theory=None, qmatoms=None, fragment=None, mm_theory=None, elif self.embedding.lower() == "mechanical" or self.embedding.lower() == "mech": self.embedding="mech" self.PC = False + elif self.embedding.lower() == "polembed_drude" or self.embedding.lower() == "drude": + self.embedding="polembed_drude" + self.PC = True else: print("Unknown embedding. Valid options are: elstat (synonyms: electrostatic, electronic), mech (synonym: mechanical)") ashexit() @@ -312,16 +315,19 @@ def __init__(self, qm_theory=None, qmatoms=None, fragment=None, mm_theory=None, if self.embedding.lower() == "elstat": print("Charges of QM atoms set to 0 (since Electrostatic Embedding):") self.ZeroQMCharges() #Modifies self.charges_qmregionzeroed - # print("length of self.charges_qmregionzeroed :", len(self.charges_qmregionzeroed)) - # TODO: make sure this works for OpenMM and for NonBondedTheory # Updating charges in MM object. self.mm_theory.update_charges(self.qmatoms,[0.0 for i in self.qmatoms]) - + elif self.embedding.lower() == "polembed_drude": + print("Polembed Drude embedding enabled.") + print("This means that QM-atoms will be zeroed for QM-MM interactions calculated by QM program") + print("But MM program will have charged defined for QM-region") + print("Not implemented yet. Exiting") + ashexit() + self.ZeroQMCharges() #Modifies self.charges_qmregionzeroed # Also removing QM-MM Coulomb interaction exceptions in OpenMM if self.mm_theory_name == "OpenMMTheory": # Deleting Coulomb exception interactions involving QM and MM atoms self.mm_theory.delete_exceptions(self.qmatoms) - elif self.embedding.lower() == "pbcmm-elstat": print("PBC Electrostatic embedding enabled.") print("This means that QM-atoms will be zeroed for QM-MM interactions calculated by QM program") @@ -775,6 +781,9 @@ def run(self, current_coords=None, elems=None, Grad=False, numcores=1, exit_afte # MM-program thus double-counts (SR QM-QM and SR QM-MM) and we need subtractive corrections return self.elstat_run(current_coords=current_coords, elems=elems, Grad=Grad, numcores=numcores, exit_after_customexternalforce_update=exit_after_customexternalforce_update, label=label, charge=charge, mult=mult) + elif self.embedding.lower() == "polembed_drude": + return self.elstat_run(current_coords=current_coords, elems=elems, Grad=Grad, numcores=numcores, exit_after_customexternalforce_update=exit_after_customexternalforce_update, + label=label, charge=charge, mult=mult) else: print("Unknown embedding. Exiting") ashexit() @@ -1108,7 +1117,7 @@ def runprep(self, current_coords): # If no linkatoms then use original self.qmelems self.current_qmelems = self.qmelems # If no linkatoms then self.pointcharges are just original charges with QM-region zeroed - if self.embedding.lower() == "elstat": + if self.embedding.lower() == "elstat" or self.embedding.lower() == "polembed_drude": self.pointcharges = [self.charges_qmregionzeroed[i] for i in self.mmatoms] # NOTE: Now we have updated MM-coordinates (if doing linkatoms, with dipolecharges etc) and updated mm-charges (more, due to dipolecharges if linkatoms) @@ -1124,7 +1133,7 @@ def runprep(self, current_coords): self.pointcharges_original=copy.copy(self.pointcharges) # Initialize QM_PC_gradient for efficiency - if self.embedding.lower() == "elstat": + if self.embedding.lower() == "elstat" or self.embedding.lower() == "polembed_drude": self.QM_PC_gradient = np.zeros((len(self.allatoms), 3)) print_time_rel(init_time_runprep, modulename='runprep', moduleindex=3, currprintlevel=self.printlevel, currthreshold=2) @@ -1291,7 +1300,7 @@ def elstat_run(self, current_coords=None, elems=None, Grad=False, numcores=1, ex self.QMenergy = QMenergy #No TruncPC approximation active. No change to original QM and PCgradient from QMcode self.QMgradient_wo_linkatoms = QMgradient_wo_linkatoms - if self.embedding.lower() == "elstat": + if self.embedding.lower() == "elstat" or self.embedding.lower() == "polembed_drude": self.PCgradient = PCgradient # Populatee QM_PC gradient (has full system size) diff --git a/ash/modules/module_coords.py b/ash/modules/module_coords.py index 4840b43be..4d4e0e99f 100644 --- a/ash/modules/module_coords.py +++ b/ash/modules/module_coords.py @@ -82,7 +82,7 @@ def calculate_reaction_energy(self): # ASH Fragment class class Fragment: - def __init__(self, fragments=None, coordsstring=None, fragfile=None, databasefile=None, xyzfile=None, pdbfile=None, grofile=None, + def __init__(self, fragments=None, coordsstring=None, fragfile=None, databasefile=None, xyzfile=None, pdbfile=None, pdbxfile=None, grofile=None, amber_inpcrdfile=None, amber_prmtopfile=None, trexiofile=None, smiles=None, chemshellfile=None, coords=None, elems=None, connectivity=None, atom=None, diatomic=None, diatomic_bondlength=None, bondlength=None, @@ -141,7 +141,6 @@ def __init__(self, fragments=None, coordsstring=None, fragfile=None, databasefil ############################## # NOW PROCESSING INPUT DATA ############################## - # Lists of elements and coordinates provided if coords is not None: # Adding coords as list of lists (or np.array). Conversion to numpy array @@ -221,6 +220,10 @@ def __init__(self, fragments=None, coordsstring=None, fragfile=None, databasefil self.label = pdbfile.split('/')[-1].split('.')[0] self.read_pdbfile_openmm(pdbfile) #self.read_pdbfile_old(pdbfile, conncalc=False, use_atomnames_as_elements=use_atomnames_as_elements) + # PDBX-file + elif pdbxfile is not None: + self.label = pdbxfile.split('/')[-1].split('.')[0] + self.read_pdbxfile(pdbxfile) # GROMACS GRO-file elif grofile is not None: self.label = grofile.split('/')[-1].split('.')[0] @@ -582,6 +585,31 @@ def read_pdbfile_openmm(self,filename): ashexit() pdb = openmm.app.PDBFile(filename) self.coords = np.array([[i.x*10,i.y*10,i.z*10] for i in pdb.positions]) + self.elems=[] + print(pdb.topology) + for atom in pdb.topology.atoms(): + + try: + self.elems.append(atom.element.symbol) + except AttributeError: + print("Warning: could not fully parse element information from PDB-topology for atom:", atom) + print("This may be a virtual site. Adding 'M' as dummy element for this atom.") + self.elems.append('M') + + # Topology + self.pdb_topology = pdb.topology + + # Reading PDBx/mmCIF file using OpenMM + def read_pdbxfile(self, filename): + if self.printlevel >= 2: + print("read_pdbxfile: Reading coordinates from PDBX file '{}' into fragment.".format(filename)) + try: + import openmm.app + except ImportError: + print("Error: OpenMM library not found. ASH requires OpenMM library to read PDB files.") + ashexit() + pdb = openmm.app.PDBxFile(filename) + self.coords = np.array([[i.x*10,i.y*10,i.z*10] for i in pdb.positions]) self.elems = [atom.element.symbol for atom in pdb.topology.atoms()] # Topology @@ -877,6 +905,49 @@ def write_pdbfile_openmm(self,filename="Fragment", calc_connectivity=False, pdb_ print(f"Wrote PDB-file: {filename}") return filename + # Write PDBX-file via OpenMM + def write_pdbxfile(self,filename="Fragment", calc_connectivity=False, pdb_topology=None, + skip_connectivity=False, resname="MOL"): + print("write_pdbxfile_openmm\n") + try: + import openmm.app + except ImportError: + print("Error: OpenMM library not found. ASH requires OpenMM library to write PDB files.") + ashexit() + + #Adding extension + if '.cif' not in filename: + filename += ".cif" + + if pdb_topology is not None: + print("Using input pdb_topology") + self.pdb_topology = pdb_topology + elif self.pdb_topology is None: + print("Warning: ASH Fragment has no PDB-file topology defined (required for PDB-file writing)") + print("Now defining new topology from scratch") + if pdb_topology is None: + self.define_topology(resname=resname) #Creates self.pdb_topology + else: + print("Using pdbtopology found in ASH fragment") + + # Before writing PDB-file, request connectivity calculation so that we get correct CONECT lines for non-biomolecules + if calc_connectivity is True: + print("Connectivity calculation requested for Fragment") + connectivity_dict = get_connected_atoms_dict(self.coords,self.elems, 1.0,0.1) + print("Adding connectivity to PDB topology") + ash.interfaces.interface_OpenMM.openmm_add_bonds_to_topology(self.pdb_topology,connectivity_dict) + + # If no_connectivity is True, we skip adding connectivity to PDB-file + if skip_connectivity is True: + print("skip_connectivity True: this will not write connectivity lines to PDB-file") + # solute_resname= list(self.pdb_topology.residues())[0].name + print("Deleting molecule bond information") + # Setting list of bonds to empty list + self.pdb_topology._bonds=[] + openmm.app.PDBxFile.writeFile(self.pdb_topology, self.coords, file=open(f"{filename}", 'w')) + print(f"Wrote PDBx-file: {filename}") + return filename + def write_xyzfile(self, xyzfilename="Fragment-xyzfile.xyz", writemode='w', write_chargemult=True, write_energy=True): with open(xyzfilename, writemode) as ofile: @@ -1124,7 +1195,13 @@ def reformat_list_to_array(l): # Can also convert atomic-number (isatomnum flag) def reformat_element(elem, isatomnum=False): if isatomnum is True: - el_correct = ash.dictionaries_lists.element_dict_atnum[elem].symbol + try: + el_correct = ash.dictionaries_lists.element_dict_atnum[elem].symbol + except KeyError: + print("Element-string: {} not found in element-dictionary!".format(elem)) + print("This is not a valid element as defined in ASH source-file: dictionaries_lists.py") + print("Fix element-information in coordinate-file.") + ashexit() else: try: el_correct = ash.dictionaries_lists.element_dict_atname[elem.lower()].symbol @@ -2482,7 +2559,11 @@ def totmasslist(ellist): for e in ellist: try: atcharge = int(elematomnumbers[e.lower()]) - atmass = atommasses[atcharge - 1] + if atcharge == 0: + print("Warning: element '{}' has atomic number 0. This is likely a dummy atom. Using mass of 0.0".format(e)) + atmass=0.0 + else: + atmass = atommasses[atcharge - 1] except KeyError: atmass = 0.0 if warning_issued is False: @@ -2501,7 +2582,11 @@ def list_of_masses(ellist): for e in ellist: try: atcharge = int(elematomnumbers[e.lower()]) - atmass = atommasses[atcharge - 1] + if atcharge == 0: + print("Warning: element '{}' has atomic number 0. This is likely a dummy atom. Using mass of 0.0".format(e)) + atmass=0.0 + else: + atmass = atommasses[atcharge - 1] except KeyError: atmass = 0.0 if warning_issued is False: @@ -2540,7 +2625,7 @@ def flexible_align_pdb(pdbfileA, pdbfileB, rotate_only=False, translate_only=Fal #Write PDBfile. PDB-info will have been read and stored fragmentA.coords=newfragA.coords #Replacing coords in original fragmentA - fragmentA.write_pdbfile(filename=f"{pdbfileA.replace('.pdb','')}_aligned") #Now write out + fragmentA.write_pdbfile_openmm(filename=f"{pdbfileA.replace('.pdb','')}_aligned") #Now write out #For ASH fragments def flexible_align(fragmentA, fragmentB, rotate_only=False, translate_only=False, reordering=False, reorder_method='brute', subset=None): diff --git a/ash/modules/module_freq.py b/ash/modules/module_freq.py index e8115a608..9c0e512a4 100644 --- a/ash/modules/module_freq.py +++ b/ash/modules/module_freq.py @@ -21,7 +21,7 @@ #TODO: IR/Raman intensities def AnFreq(fragment=None, theory=None, charge=None, mult=None, temp=298.15, masses=None, pressure=1.0, QRRHO=True, QRRHO_method='Grimme', QRRHO_omega_0=100, printlevel=2, - scaling_factor=1.0, symmetry_number=None): + scaling_factor=1.0, symmetry_number=None, rotmode_threshold=1e-4): module_init_time=time.time() print(BC.WARNING, BC.BOLD, "------------ANALYTICAL FREQUENCIES-------------", BC.END) @@ -30,7 +30,7 @@ def AnFreq(fragment=None, theory=None, charge=None, mult=None, temp=298.15, mass ashexit() # Checking for linearity. Determines how many Trans+Rot modes - if detect_linear(coords=fragment.coords,elems=fragment.elems) is True: + if detect_linear(coords=fragment.coords,elems=fragment.elems, threshold=rotmode_threshold) is True: TRmodenum=5 else: TRmodenum=6 @@ -64,7 +64,7 @@ def AnFreq(fragment=None, theory=None, charge=None, mult=None, temp=298.15, mass # Would ensure completely correct masses at least # For now grabbing directly from theory object # Tested with pyscf, ORCA - + IR_intens_values=None try: IR_intens_values=theory.ir_intensities if len(IR_intens_values) == 0: @@ -86,7 +86,7 @@ def AnFreq(fragment=None, theory=None, charge=None, mult=None, temp=298.15, mass print("Normal mode composition factors by element") printfreqs_and_nm_elem_comps(frequencies,fragment,evectors,hessatoms=hessatoms,TRmodenum=TRmodenum) thermodict = thermochemcalc(frequencies,hessatoms, fragment, mult, temp=temp,pressure=pressure, QRRHO=QRRHO, QRRHO_omega_0=QRRHO_omega_0, - symmetry_number=symmetry_number) + symmetry_number=symmetry_number, rotmode_threshold=rotmode_threshold) # Add Hessian to fragment and write to file fragment.hessian=hessian @@ -113,7 +113,8 @@ def AnFreq(fragment=None, theory=None, charge=None, mult=None, temp=298.15, mass # Numerical frequencies function # ORCA uses 0.005 Bohr = 0.0026458861 Ang, CHemshell uses 0.01 Bohr = 0.00529 Ang def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displacement=0.005, hessatoms=None, numcores=1, runmode='serial', - temp=298.15, pressure=1.0, hessatoms_masses=None, printlevel=1, QRRHO=True, QRRHO_method='Grimme', QRRHO_omega_0=100, Raman=False, + temp=298.15, pressure=1.0, hessatoms_masses=None, printlevel=1, QRRHO=True, QRRHO_method='Grimme', QRRHO_omega_0=100, + IR=True, Raman=False, rotmode_threshold=1e-4, scaling_factor=1.0, symmetry_number=None, force_projection=None): module_init_time=time.time() print(BC.WARNING, BC.BOLD, "------------NUMERICAL FREQUENCIES-------------", BC.END) @@ -181,7 +182,7 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa print("Check input masses!",BC.END) ashexit() # Checking for linearity. Determines how many Trans+Rot modes - if detect_linear(coords=fragment.coords,elems=fragment.elems) is True: + if detect_linear(coords=fragment.coords,elems=fragment.elems, threshold=rotmode_threshold) is True: TRmodenum=5 else: TRmodenum=6 @@ -335,11 +336,13 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa displacement_grad_dictionary[stringlabel] = gradient #Grabbing dipole moment if available - try: - displacement_dm = theory.get_dipole_moment() - displacement_dipole_dictionary[stringlabel] = displacement_dm - except: - pass + if IR is True: + try: + displacement_dm = theory.get_dipole_moment() + displacement_dipole_dictionary[stringlabel] = displacement_dm + except: + pass + #Grabbing polarizability tensor if requested if Raman is True: try: @@ -411,9 +414,10 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa #original_grad=get_partial_matrix(allatoms, hessatoms, displacement_grad_dictionary['Originalgeo']) original_grad_1d = np.ravel(original_grad) #IR intensities if dipoles available - if len(displacement_dipole_dictionary) > 0: - original_dipole = np.array(displacement_dipole_dictionary['Originalgeo']) - #print("original_dipole:",original_dipole) + if IR is True: + if len(displacement_dipole_dictionary) > 0: + original_dipole = np.array(displacement_dipole_dictionary['Originalgeo']) + #print("original_dipole:",original_dipole) #Raman if requested if Raman is True: if len(displacement_polarizability_dictionary) > 0: @@ -438,15 +442,16 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa hessian[hessindex,:]=Hessrow grad_pos_1d=0 #IR intensities if dipoles available - if len(displacement_dipole_dictionary) > 0: - # Make sure it's not a dict of None's - if any(value is None for value in displacement_dipole_dictionary.values()): - #print("None values in displacement_dipole_dictionary. Skipping IR") - pass - elif len(displacement_dipole_dictionary[lookup_string_pos]) > 0: - disp_dipole = np.array(displacement_dipole_dictionary[lookup_string_pos]) - dd_deriv = (disp_dipole - original_dipole)/displacement_bohr - dipole_derivs[hessindex,:] = dd_deriv + if IR is True: + if len(displacement_dipole_dictionary) > 0: + # Make sure it's not a dict of None's + if any(value is None for value in displacement_dipole_dictionary.values()): + #print("None values in displacement_dipole_dictionary. Skipping IR") + pass + elif len(displacement_dipole_dictionary[lookup_string_pos]) > 0: + disp_dipole = np.array(displacement_dipole_dictionary[lookup_string_pos]) + dd_deriv = (disp_dipole - original_dipole)/displacement_bohr + dipole_derivs[hessindex,:] = dd_deriv #Raman if requested if Raman is True: if len(displacement_polarizability_dictionary) > 0: @@ -486,7 +491,7 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa grad_neg_1d=0 #IR intensities if dipoles available - if len(displacement_dipole_dictionary) > 0: + if IR is True and len(displacement_dipole_dictionary) > 0: # Make sure it's not a dict of None's if any(value is None for value in displacement_dipole_dictionary.values()): #print("None values in displacement_dipole_dictionary. Skipping IR") @@ -537,17 +542,18 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa #Evectors: eigenvectors of the mass-weighed Hessian #Normal modes: unweighted - frequencies, nmodes, evectors, mode_order = diagonalizeHessian(hesscoords,hessian,hessmasses,hesselems,TRmodenum=TRmodenum,projection=projection) + frequencies, nmodes, evectors, mode_order = diagonalizeHessian(hesscoords,hessian,hessmasses,hesselems, + TRmodenum=TRmodenum,projection=projection, rotmode_threshold=rotmode_threshold) print("Diagonalization of frequencies complete") print("Now scaling frequencies by scaling factor:", scaling_factor) frequencies = scaling_factor * np.array(frequencies) # IR intensities if dipoles available - if np.any(dipole_derivs): - dipole_derivs = dipole_derivs[mode_order] - IR_intens_values = calc_IR_Intensities(hessmasses,evectors,dipole_derivs) - else: - IR_intens_values = None + IR_intens_values = None + if IR is True: + if np.any(dipole_derivs): + dipole_derivs = dipole_derivs[mode_order] + IR_intens_values = calc_IR_Intensities(hessmasses,evectors,dipole_derivs) # Raman activities if polarizabilities available if Raman is True: @@ -577,7 +583,7 @@ def NumFreq(fragment=None, theory=None, charge=None, mult=None, npoint=2, displa # Get and print out thermochemistry thermodict = thermochemcalc(frequencies,hessatoms, fragment, mult, temp=temp,pressure=pressure, QRRHO=QRRHO, QRRHO_method=QRRHO_method, QRRHO_omega_0=QRRHO_omega_0, - symmetry_number=symmetry_number) + symmetry_number=symmetry_number, rotmode_threshold=rotmode_threshold) # Write Hessian to file write_hessian(hessian,hessfile="Hessian") @@ -624,9 +630,9 @@ def get_partial_matrix(matrix,hessatoms): # Diagonalize Hessian from input Hessian, masses and element-strings -def diagonalizeHessian(coords,hessian, masses, elems, projection=True, TRmodenum=None, LargeImagFreqThreshold=-100): +def diagonalizeHessian(coords,hessian, masses, elems, projection=True, TRmodenum=None, + LargeImagFreqThreshold=-100, rotmode_threshold=1e-4): print("\nDiagonalizing Hessian") - numatoms=len(elems) atomlist = [] for i, j in enumerate(elems): atomlist.append(str(j) + '-' + str(i)) @@ -634,8 +640,7 @@ def diagonalizeHessian(coords,hessian, masses, elems, projection=True, TRmodenum # Projecting out translations and rotations if projection is True: print("Projection of out rotational and translational modes active!") - vfreqs,evectors,nmodes = project_rot_and_trans(coords,masses,hessian) - + vfreqs,evectors,nmodes = project_rot_and_trans(coords,masses,hessian,rotmode_threshold=rotmode_threshold) # Adding TRmodes zeros to vfreqs list for i in range(0,TRmodenum): vfreqs = np.insert(vfreqs,0,0.0) @@ -736,7 +741,7 @@ def printfreqs(vfreq,numatoms,TRmodenum=6, intensities=None, Raman_activities=No print("VIBRATIONAL FREQUENCY SUMMARY") print("-"*40) if intensities is None: - print("No IR intensities were calculated (dipoles not available in QM-program interface). Setting values to 0.0.") + print("No IR intensities were calculated. Setting values to 0.0.") if Raman_activities is None: print("No Raman activities were calculated (polarizabilities not available in QM-program interface). Setting values to 0.0.") print("Note: imaginary modes shown as negative") @@ -806,7 +811,7 @@ def old_printfreqs(vfreq,numatoms,TRmodenum=6): # def thermochemcalc(vfreq,atoms,fragment, multiplicity, temp=298.15,pressure=1.0, QRRHO=True, QRRHO_method='Grimme', QRRHO_omega_0=100, - use_full_geo_in_rotational_analysis=True, symmetry_number=None): + use_full_geo_in_rotational_analysis=True, symmetry_number=None, rotmode_threshold=1e-4): module_init_time=time.time() """[summary] @@ -833,7 +838,7 @@ def thermochemcalc(vfreq,atoms,fragment, multiplicity, temp=298.15,pressure=1.0, TRmodenum=5 else: print("System size > 2, checking if linear") - linearcheck = detect_linear(fragment) + linearcheck = detect_linear(fragment, threshold=rotmode_threshold) if linearcheck is True: print("Structure is linear. 5 translational+rotational modes present") moltype="linear" @@ -1070,7 +1075,7 @@ def thermochemcalc(vfreq,atoms,fragment, multiplicity, temp=298.15,pressure=1.0, #From Hess-tool.py: Copied 13 May 2020 #Print dummy ORCA outputfile using coordinates and normal modes. Used for visualization of modes in Chemcraft -def printdummyORCAfile(elems,coords,vfreq,evectors,nmodes,hessfile): +def printdummyORCAfile(elems,coords,vfreq,evectors,nmodes,hessfile, rotmode_threshold=1e-4): orca_header = """ ***************** * O R C A * ***************** @@ -1088,7 +1093,7 @@ def printdummyORCAfile(elems,coords,vfreq,evectors,nmodes,hessfile): CARTESIAN COORDINATES (ANGSTROEM) ---------------------------------""" #Checking for linearity here. - if detect_linear(coords=coords,elems=elems) == True: + if detect_linear(coords=coords,elems=elems, threshold=rotmode_threshold) is True: TRmodenum=5 else: TRmodenum=6 @@ -2188,7 +2193,7 @@ def clean_frequencies(freqs): -def project_rot_and_trans(coords,mass,Hessian): +def project_rot_and_trans(coords,mass,Hessian, rotmode_threshold=1e-4): mass = np.array(mass) coords = np.array(coords)*ash.constants.ang2bohr coords = coords.copy().reshape(-1, 3) @@ -2218,9 +2223,11 @@ def project_rot_and_trans(coords,mass,Hessian): RotDOF = 0 for i in range(3): print("Ivals[i]:", Ivals[i]) - if abs(Ivals[i]) > 1.0e-5: + if abs(Ivals[i]) > rotmode_threshold: RotDOF += 1 TR_DOF = 3 + RotDOF + print("TR_DOF:", TR_DOF) + #exit() if TR_DOF not in (5, 6): print("Unexpected number of trans+rot DOF: {TR_DOF} not in (5, 6)") diff --git a/ash/modules/module_machine_learning.py b/ash/modules/module_machine_learning.py index 3c675c38d..fcab6c1b4 100644 --- a/ash/modules/module_machine_learning.py +++ b/ash/modules/module_machine_learning.py @@ -14,6 +14,45 @@ # Collection of functions related to machine learning and data analysis # Also helper tools for Torch and MLatom interfaces +# Helper that writes a training-data file in the MACE extended-XYZ format. +# atom_energies/energies are expected in Hartree and gradients in Hartree/Bohr; +# all are converted to MACE units (eV and eV/Å) on the fly. +def write_mace_xyz_file(filename, atom_energies, energies, gradients, fragments, Grad=True): + #TODO: Nmols, comp, molindex ? + Nmols="1" + comp="xxx" + molindex=0 + with open(filename, "w") as mace_file: + print(f"Writing isolated atom reference energies to {filename} ....") + for el, an_at in atom_energies.items(): + en_ev = an_at * 27.211386245988 + mace_file.write("1\n") + if Grad: + mace_file.write(f"Properties=species:S:1:pos:R:3:forces_REF:R:3 config_type=IsolatedAtom energy_REF={en_ev} pbc='F F F'\n") + mace_file.write(f"{el:2s}{0.0:17.8f}{0.0:17.8f}{0.0:17.8f}" + f"{-0.0:17.8f}{-0.0:17.8f}{-0.0:17.8f}\n") + else: + mace_file.write(f"Properties=species:S:1:pos:R:3 config_type=IsolatedAtom energy_REF={en_ev} pbc='F F F'\n") + mace_file.write(f"{el:2s}{0.0:17.8f}{0.0:17.8f}{0.0:17.8f}\n") + + for i in range(len(energies)): + # Converting energy to eV + frag = fragments[i] + energy_ev = energies[i]*27.211386245988 + mace_file.write(f"{frag.numatoms}\n") + if Grad: + force = -1 * np.array(gradients[i]) * 51.42206747 + mace_file.write(f"Properties=species:S:1:pos:R:3:molID:I:1:forces_REF:R:3 Nmols={Nmols} Comp={comp} energy_REF={energy_ev} pbc='F F F'\n") + for j in range(frag.numatoms): # Bug fix: inner loop variable was i, now j + mace_file.write(f"{frag.elems[j]:2s}{frag.coords[j][0]:17.8f}{frag.coords[j][1]:17.8f}{frag.coords[j][2]:17.8f}" + f"{molindex:9d}{force[j][0]:17.8f}{force[j][1]:17.8f}{force[j][2]:17.8f}\n") + else: + # Write without forces + mace_file.write(f"Properties=species:S:1:pos:R:3:molID:I:1 Nmols={Nmols} Comp={comp} energy_REF={energy_ev} pbc='F F F'\n") + for j in range(frag.numatoms): # Bug fix: inner loop variable was i, now j + mace_file.write(f"{frag.elems[j]:2s}{frag.coords[j][0]:17.8f}{frag.coords[j][1]:17.8f}{frag.coords[j][2]:17.8f}" + f"{molindex:9d}\n") + # Function to create ML training data given XYZ-files and 2 ASH theories def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=None, xyz_files=None, num_snapshots=None, random_snapshots=True, dcd_pdb_topology=None, nth_frame_in_traj=1, printlevel=2, @@ -194,7 +233,10 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No os.chdir('..') # Remove old files if present - for f in ["train_data.xyz", "train_data.energies", "train_data.gradients", "train_data_mace.xyz"]: + for f in ["train_data.xyz", "train_data.energies", "train_data.gradients", "train_data_mace.xyz", + "train_data_theory1.energies", "train_data_theory2.energies", + "train_data_theory1.gradients", "train_data_theory2.gradients", + "train_data_mace_theory1.xyz", "train_data_mace_theory2.xyz"]: try: os.remove(f) except: @@ -205,10 +247,16 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No gradients=[] fragments=[] labels=[] + # When delta-learning (two theories), also keep the individual theory contributions + energies_theory1=[] + energies_theory2=[] + gradients_theory1=[] + gradients_theory2=[] - # Removing + # Removing theory_1.cleanup() - theory_2.cleanup() + if theory_2 is not None: + theory_2.cleanup() if runmode=="serial": print("Runmode is serial!") @@ -243,6 +291,12 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No energy = result_2.energy - result_1.energy if Grad is True: gradient = result_2.gradient - result_1.gradient + # Keep individual theory contributions for separate output files + energies_theory1.append(result_1.energy) + energies_theory2.append(result_2.energy) + if Grad is True: + gradients_theory1.append(result_1.gradient) + gradients_theory2.append(result_2.gradient) else: energy = result_1.energy if Grad is True: @@ -291,6 +345,9 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No if delta is True: energy = results_theory2.energies_dict[l] - results_theory1.energies_dict[l] print("energy:", energy) + # Keep individual theory contributions for separate output files + energies_theory1.append(results_theory1.energies_dict[l]) + energies_theory2.append(results_theory2.energies_dict[l]) else: energy = results_theory1.energies_dict[l] @@ -302,11 +359,16 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No if Grad: if delta is True: gradient = results_theory2.gradients_dict[l] - results_theory1.gradients_dict[l] + gradients_theory1.append(results_theory1.gradients_dict[l]) + gradients_theory2.append(results_theory2.gradients_dict[l]) else: gradient = results_theory1.gradients_dict[l] gradients.append(gradient) # Calculate energies for atoms + # Per-theory isolated-atom reference energies, needed to write per-theory MACE files + energies_atoms_dict_theory1={} + energies_atoms_dict_theory2={} if energies_atoms_dict is None: print("\nNow calculating isolated atom reference energies for each element in the training set") energies_atoms_dict={} @@ -322,6 +384,7 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No theory_1.cleanup() result_1 = Singlepoint(theory=theory_1, fragment=atomfrag, printlevel=0, result_write_to_disk=False) + energies_atoms_dict_theory1[uniq_el] = result_1.energy if delta is True: theory_2.printlevel=0 # Running theory 2 @@ -329,6 +392,7 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No theory_2.cleanup() result_2 = Singlepoint(theory=theory_2, fragment=atomfrag, printlevel=0, result_write_to_disk=False) + energies_atoms_dict_theory2[uniq_el] = result_2.energy # Delta energy atomenergy = result_2.energy - result_1.energy else: @@ -364,52 +428,62 @@ def create_ML_training_data(xyz_dir=None, dcd_trajectory=None, xyz_trajectory=No gradients_file.write(f"{g[0]:10.7f} {g[1]:10.7f} {g[2]:10.7f}\n") gradients_file.close() + # When delta-learning (two theories), also write the individual theory + # energies and gradients to separate files alongside the delta data above. + if delta is True: + # Per-theory energy files + for theoryname, theory_energies in [("theory1", energies_theory1), + ("theory2", energies_theory2)]: + with open(f"train_data_{theoryname}.energies", "w") as f_energies: + for energy in theory_energies: + f_energies.write(f"{energy}\n") + + # Per-theory gradient files + if Grad: + for theoryname, theory_gradients in [("theory1", gradients_theory1), + ("theory2", gradients_theory2)]: + with open(f"train_data_{theoryname}.gradients", "w") as f_gradients: + for frag, grad in zip(fragments, theory_gradients): + f_gradients.write(f"{frag.numatoms}\n") + f_gradients.write(f"gradient {frag.label} \n") + for g in grad: + f_gradients.write(f"{g[0]:10.7f} {g[1]:10.7f} {g[2]:10.7f}\n") + print("\nNow writing data in MACE-format with energies in units of eV and forces in eV/Å") print("Fragments labels:",[frag.label for frag in fragments]) print("energies:", energies) - # Write data file that MACE uses - - with open("train_data_mace.xyz", "w") as mace_file: - print("Writing isolated atom reference energies....") - for el, an_at in energies_atoms_dict.items(): - en_ev = an_at * 27.211386245988 - mace_file.write("1\n") - if Grad: - mace_file.write(f"Properties=species:S:1:pos:R:3:forces_REF:R:3 config_type=IsolatedAtom energy_REF={en_ev} pbc='F F F'\n") - mace_file.write(f"{el:2s}{0.0:17.8f}{0.0:17.8f}{0.0:17.8f}" - f"{-0.0:17.8f}{-0.0:17.8f}{-0.0:17.8f}\n") - else: - mace_file.write(f"Properties=species:S:1:pos:R:3 config_type=IsolatedAtom energy_REF={en_ev} pbc='F F F'\n") - mace_file.write(f"{el:2s}{0.0:17.8f}{0.0:17.8f}{0.0:17.8f}\n") - #TODO: Nmols, comp, molindex ? - Nmols="1" - comp="xxx" - molindex=0 - - for i in range(len(energies)): - # Converting energy to eV - frag = fragments[i] - energy_ev = energies[i]*27.211386245988 - mace_file.write(f"{frag.numatoms}\n") - if Grad: - force = -1 * np.array(gradients[i]) * 51.42206747 - mace_file.write(f"Properties=species:S:1:pos:R:3:molID:I:1:forces_REF:R:3 Nmols={Nmols} Comp={comp} energy_REF={energy_ev} pbc='F F F'\n") - for j in range(frag.numatoms): # Bug fix: inner loop variable was i, now j - mace_file.write(f"{frag.elems[j]:2s}{frag.coords[j][0]:17.8f}{frag.coords[j][1]:17.8f}{frag.coords[j][2]:17.8f}" - f"{molindex:9d}{force[j][0]:17.8f}{force[j][1]:17.8f}{force[j][2]:17.8f}\n") - else: - # Write without forces - mace_file.write(f"Properties=species:S:1:pos:R:3:molID:I:1 Nmols={Nmols} Comp={comp} energy_REF={energy_ev} pbc='F F F'\n") - for j in range(frag.numatoms): # Bug fix: inner loop variable was i, now j - mace_file.write(f"{frag.elems[j]:2s}{frag.coords[j][0]:17.8f}{frag.coords[j][1]:17.8f}{frag.coords[j][2]:17.8f}" - f"{molindex:9d}\n") + # Write combined data file that MACE uses (delta data when two theories are used, + # otherwise the single-theory data) + write_mace_xyz_file("train_data_mace.xyz", energies_atoms_dict, energies, gradients, + fragments, Grad=Grad) + + # When delta-learning (two theories), also write a MACE-format file for each + # individual theory level, using that theory's own atomic reference energies. + if delta is True: + if energies_atoms_dict_theory1 and energies_atoms_dict_theory2: + write_mace_xyz_file("train_data_mace_theory1.xyz", energies_atoms_dict_theory1, + energies_theory1, gradients_theory1, fragments, Grad=Grad) + write_mace_xyz_file("train_data_mace_theory2.xyz", energies_atoms_dict_theory2, + energies_theory2, gradients_theory2, fragments, Grad=Grad) + else: + print("Warning: per-theory isolated-atom reference energies are not available " + "(user-provided energies_atoms_dict only holds the combined values).") + print("Skipping per-theory MACE files train_data_mace_theory1.xyz / train_data_mace_theory2.xyz") print("All done! Files created:\ntrain_data.xyz\ntrain_data.energies\ntrain_data_mace.xyz") if Grad: print("train_data.gradients") + if delta is True: + print("train_data_theory1.energies\ntrain_data_theory2.energies") + if Grad: + print("train_data_theory1.gradients\ntrain_data_theory2.gradients") + if energies_atoms_dict_theory1 and energies_atoms_dict_theory2: + print("train_data_mace_theory1.xyz\ntrain_data_mace_theory2.xyz") print("Number of user-chosen snapshots:", num_snapshots) print("Number of successfully generated datapoints:", len(energies)) + return len(energies) + # Print statistics for dict with statistics for many models # Assumes a dictionary with modelfilenames as keys and statistics_dict_forDB as values @@ -741,3 +815,60 @@ def move_chosen_files(chosen,dirname): # WORKFLOW FUNCTIONS FOR TRAINING ####################################### +# Distance matrix +def distance_matrix(fragment): + from scipy.spatial.distance import cdist + return cdist(fragment.coords, fragment.coords) + +# Geometric bond-order matrix +def geomBO(fragment, alpha=6.0, cutoff_factor=1.5): + + from scipy.spatial.distance import cdist + # Distance matrix + D = cdist(fragment.coords, fragment.coords) + + # ============================================================ + # Geometry-based bond-order matrix + # B_ij = exp[-alpha * (D_ij / Rij - 1)] + # Rij = covalent radius sum + # ============================================================ + + N = fragment.numatoms + B = np.zeros((N, N)) + from ash.modules.module_coords import eldict_covrad + for i in range(N): + for j in range(i + 1, N): + Ri = eldict_covrad[fragment.elems[i]] + Rj = eldict_covrad[fragment.elems[j]] + Rij = Ri + Rj + dij = D[i, j] + # optional cutoff + if dij > cutoff_factor * Rij: + value = 0.0 + else: + value = np.exp( + -alpha * (dij / Rij - 1.0) + ) + B[i, j] = value + B[j, i] = value + return B + +# Coulomb matrix +def coulomb_matrix(fragment): + + from scipy.spatial.distance import cdist + + # Distance matrix + D = cdist(fragment.coords, fragment.coords) + N = fragment.numatoms + C = np.zeros((N, N)) + Z = fragment.nuc_charges + for i in range(N): + for j in range(N): + if i == j: + # diagonal term + C[i, i] = 0.5 * Z[i] ** 2.4 + else: + # avoid divide-by-zero + C[i, j] = Z[i] * Z[j] / D[i, j] + return C diff --git a/ash/modules/module_singlepoint.py b/ash/modules/module_singlepoint.py index 76f2e5b99..a7debfeb6 100644 --- a/ash/modules/module_singlepoint.py +++ b/ash/modules/module_singlepoint.py @@ -224,7 +224,7 @@ def Singlepoint_fragments(theory=None, fragments=None, stoichiometry=None, relat #Printing reaction energy if stoichiometry was provided if stoichiometry != None: print("Stoichiometry provided:", stoichiometry) - r = ReactionEnergy(list_of_energies=energies, stoichiometry=stoichiometry, list_of_fragments=fragments, unit='kcal/mol', label='ΔE') + r = ReactionEnergy(list_of_energies=energies, stoichiometry=stoichiometry, list_of_fragments=fragments, unit=unit, label='ΔE') result.reaction_energy = r[0] result.write_to_disk(filename="ASH_SP_fragments.result") print_time_rel(module_init_time, modulename='Singlepoint_fragments', moduleindex=1) diff --git a/pyproject.toml b/pyproject.toml index ea046e213..6379c3852 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "ash" -version = "0.95" +version = "0.99" description = "ASH is a multiscale, multitheory modelling program" authors = [ {name = "R. Bjornsson", email = "ragnar.bjornsson@gmail.com" }