Merge pull request #505 from RagnarB83/NEW

Bugfixes and some small new features

Photoelectron: CCSD(T) option, OODFT moread and outputfiles fix
DLFIND interface: bugfix for TSOpt
create_ML_training_data: more data files written, some restructuring (Madhur)
PDBx/mmCIF file support for Fragment and OpenMMTHeory
flexible_align_pdb: fix
fixes for elements dictionaries: lanthanides and actinides, also dummy atom support for certain FFs
pySCF: skala use fix
kspies: minor fixes
crest interface: support for socket
CFOUR: bugfix for NEB runs (Madhur)
SellaIRC option (undocumented)

Author: RagnarB83

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 <https://onlinelibrary.wiley.com/doi/10.1002/jcc.70359>`_
-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 <https://onlinelibrary.wiley.com/doi/10.1002/jcc.70359>`_
+R. Bjornsson, J. Comput. Chem 2026, 47, e70359.
 
 **In case of problems:**
 Please open an issue on Github and we will try to fix any problems as soon as possible.

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

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,

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
 

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