diff --git a/examples/job/elstat_decomp.py b/examples/job/elstat_decomp.py new file mode 100644 index 00000000..9935592c --- /dev/null +++ b/examples/job/elstat_decomp.py @@ -0,0 +1,62 @@ +from tcutility import job, results +from scm import plams + + +def load(xyz): + mol = plams.Molecule() + for line in xyz.strip().splitlines(): + el, x, y, z = line.split()[:4] + mol.add_atom(plams.Atom(symbol=el, coords=(x, y, z))) + return mol + +mol = load(""" +N 0.00000000 0.00000000 0.81411087 +H 0.47592788 -0.82433127 1.14462651 +H 0.47592788 0.82433127 1.14462651 +H -0.95185576 0.00000000 1.14462651 +B 0.00000000 0.00000000 -0.83480980 +H -0.58140445 -1.00702205 -1.13772686 +H -0.58140445 1.00702205 -1.13772686 +H 1.16280890 0.00000000 -1.13772686 +""") + +# optimize mol +with job.ADFJob() as opt_job: + opt_job.rundir = 'Elstat' + opt_job.name = 'Opt' + opt_job.molecule(mol) + opt_job.optimization() + opt_job.functional('BP86') + opt_job.basis_set('TZP') + opt_job.quality('Good') + + +with job.ADFFragmentJob(decompose_elstat=True) as frag_job: + frag_job.molecule(opt_job.output_mol_path) + frag_job.rundir = 'Elstat' + frag_job.name = 'EDA' + frag_job.functional('BP86') + frag_job.basis_set('TZP') + frag_job.quality('Good') + + frag_job.add_fragment([1, 2, 3, 4], 'Donor') + frag_job.add_fragment([5, 6, 7, 8], 'Acceptor') + +res_main = results.read('Elstat/EDA/complex_STOFIT') +print(res_main.properties.energy.elstat) + + +res_acceptor = results.read('Elstat/EDA/complex_Acceptor_NoElectrons') +print(res_acceptor.properties.energy.elstat) + + +res_donor = results.read('Elstat/EDA/complex_Donor_NoElectrons') +print(res_donor.properties.energy.elstat) + + +print('Electrostatic Potential Terms:') +print(f' e–e repulsion: {res_main.properties.energy.elstat.Vee: 10.0f} kcal/mol') +print(f' N–N repulsion: {res_main.properties.energy.elstat.Vnn: 10.0f} kcal/mol') +print(f' e(Donor)–N(Acceptor) attraction: {res_acceptor.properties.energy.elstat.Ven: 10.0f} kcal/mol') +print(f' N(Donor)–e(Acceptor) attraction: {res_donor.properties.energy.elstat.Ven: 10.0f} kcal/mol') +print(f' total: {res_main.properties.energy.elstat.total: 10.1f} kcal/mol') diff --git a/src/tcutility/job/adf.py b/src/tcutility/job/adf.py index 21bb175b..86c74353 100644 --- a/src/tcutility/job/adf.py +++ b/src/tcutility/job/adf.py @@ -263,7 +263,9 @@ def _check_job(self): class ADFFragmentJob(ADFJob): def __init__(self, *args, **kwargs): - self.childjobs = {} + self.decompose_elstat = kwargs.pop('decompose_elstat', False) + self.child_jobs = {} + # self.childjobs_no_electrons = {} super().__init__(*args, **kwargs) self.name = 'EDA' @@ -297,25 +299,25 @@ def add_fragment(self, mol: plams.Molecule, name: str = None, charge: int = 0, s # check if the atoms in the new fragment are already present in the other fragments. # if it is we should raise an error - for child in self.childjobs.values(): + for child in self.child_jobs.values(): if any((atom.symbol, atom.coords) == (myatom.symbol, myatom.coords) for atom in child._molecule for myatom in mol): log.error('An atom is present in multiple fragments.') return - name = name or f'fragment{len(self.childjobs) + 1}' - self.childjobs[name] = ADFJob(test_mode=self.test_mode) - self.childjobs[name].molecule(mol) - self.childjobs[name].charge(charge) - self.childjobs[name].spin_polarization(spin_polarization) - setattr(self, name, self.childjobs[name]) + name = name or f'fragment{len(self.child_jobs) + 1}' + self.child_jobs[name] = ADFJob(test_mode=self.test_mode) + self.child_jobs[name].molecule(mol) + self.child_jobs[name].charge(charge) + self.child_jobs[name].spin_polarization(spin_polarization) + setattr(self, name, self.child_jobs[name]) if not add_frag_to_mol: return if self._molecule is None: - self._molecule = self.childjobs[name]._molecule.copy() + self._molecule = self.child_jobs[name]._molecule.copy() else: - for atom in self.childjobs[name]._molecule.copy(): + for atom in self.child_jobs[name]._molecule.copy(): if any((atom.symbol, atom.coords) == (myatom.symbol, myatom.coords) for myatom in self._molecule): continue self._molecule.add_atom(atom) @@ -368,10 +370,10 @@ def remove_virtuals(self, frag=None, subspecies=None, nremove=None): # sum up the number of virtuals per atom in the fragment nremove = 0 - for atom in self.childjobs[frag]._molecule: + for atom in self.child_jobs[frag]._molecule: nremove += data.basis_sets.number_of_virtuals(atom.symbol, self._basis_set) # positive charge adds a virtual and negative removes a virtual - nremove += self.childjobs[frag].settings.input.ams.System.charge or 0 + nremove += self.child_jobs[frag].settings.input.ams.System.charge or 0 self.settings.input.adf.setdefault('RemoveFragOrbitals', '') self.settings.input.adf.RemoveFragOrbitals += f""" @@ -387,19 +389,19 @@ def run(self): It will also submit a calculation with the SCF iterations set to 0 in order to facilitate investigation of the field effects using PyOrb. ''' # check if the user defined fragments for this job - if len(self.childjobs) == 0: + if len(self.child_jobs) == 0: log.warn('Fragments were not specified yet, trying to read them from the xyz file ...') # if they did not define the fragments, try to guess them using the xyz-file if not self.guess_fragments(): log.error('Please specify the fragments using ADFFragmentJob.add_fragment or specify them in the xyz file.') raise - mol_str = " + ".join([formula.molecule(child._molecule) for child in self.childjobs.values()]) + mol_str = " + ".join([formula.molecule(child._molecule) for child in self.child_jobs.values()]) log.flow(f'ADFFragmentJob [{mol_str}]', ['start']) # obtain some system wide properties of the molecules - charge = sum([child.settings.input.ams.System.charge or 0 for child in self.childjobs.values()]) - unrestricted = any([(child.settings.input.adf.Unrestricted or 'no').lower() == 'yes' for child in self.childjobs.values()]) - spinpol = sum([child.settings.input.adf.SpinPolarization or 0 for child in self.childjobs.values()]) + charge = sum([child.settings.input.ams.System.charge for child_name, child in self.child_jobs.items() if not child_name.endswith('_NoElectrons')]) + unrestricted = any([(child.settings.input.adf.Unrestricted or 'no').lower() == 'yes' for child_name, child in self.child_jobs.items() if not child_name.endswith('_NoElectrons')]) + spinpol = sum([child.settings.input.adf.SpinPolarization for child_name, child in self.child_jobs.items() if not child_name.endswith('_NoElectrons')]) log.flow(f'Level: {self._functional}/{self._basis_set}') log.flow(f'Solvent: {self._solvent}') log.flow(f'Charge: {charge}', ['straight']) @@ -407,24 +409,23 @@ def run(self): log.flow(f'Spin-Polarization: {spinpol}', ['straight']) log.flow() # this job and all its children should have the same value for unrestricted - [child.unrestricted(unrestricted) for child in self.childjobs.values()] + [child.unrestricted(unrestricted) for child in self.child_jobs.values()] - # propagate the post- and preambles to the childjobs - [child.add_preamble(preamble) for preamble in self._preambles for child in self.childjobs.values()] - [child.add_postamble(postamble) for postamble in self._postambles for child in self.childjobs.values()] + # propagate the post- and preambles to the child_jobs + [child.add_preamble(preamble) for preamble in self._preambles for child in self.child_jobs.values()] + [child.add_postamble(postamble) for postamble in self._postambles for child in self.child_jobs.values()] # we now update the child settings with the parent settings # this is because we have to propagate settings such as functionals, basis sets etc. sett = self.settings.as_plams_settings() # first create a plams settings object # same for the children - child_setts = {name: child.settings.as_plams_settings() for name, child in self.childjobs.items()} + child_setts = {name: child.settings.as_plams_settings() for name, child in self.child_jobs.items()} # update the children using the parent settings [child_sett.update(sett) for child_sett in child_setts.values()] [child_sett.input.adf.pop('RemoveFragOrbitals', None) for child_sett in child_setts.values()] [child_sett.input.adf.pop('RemoveAllFragVirtuals', None) for child_sett in child_setts.values()] # same for sbatch settings - [child.sbatch(**self._sbatch) for child in self.childjobs.values()] - + [child.sbatch(**self._sbatch) for child in self.child_jobs.values()] # now set the charge, spinpol, unrestricted for the parent self.charge(charge) self.spin_polarization(spinpol) @@ -432,50 +433,101 @@ def run(self): if unrestricted: self.settings.input.adf.UnrestrictedFragments = 'Yes' + elstat_jobs = {} + # now we are going to run each child job - for i, (childname, child) in enumerate(self.childjobs.items(), start=1): - log.flow(f'Child job ({i}/{len(self.childjobs)}) {childname} [{formula.molecule(child._molecule)}]', ['split']) - log.flow(f'Charge: {child.settings.input.ams.System.charge or 0}', ['straight', 'straight']) - log.flow(f'Spin-Polarization: {child.settings.input.adf.SpinPolarization or 0}', ['straight', 'straight']) + for i, (child_name, child) in enumerate(self.child_jobs.items(), start=1): # the child name will be prepended with SP showing that it is the singlepoint calculation - child.name = f'frag_{childname}' + child.name = f'frag_{child_name}' child.rundir = j(self.rundir, self.name) - # add the path to the child adf.rkf file as a dependency to the parent job - self.settings.input.adf.fragments[childname] = j(child.workdir, 'adf.rkf') + # # add the path to the child adf.rkf file as a dependency to the parent job + self.settings.input.adf.fragments[child_name] = j(child.workdir, 'adf.rkf') + # recast the plams.Settings object into a Result object as that is what run expects + child.settings = results.Result(child_setts[child_name]) + log.flow(f'Fragment ({i}/{len(self.child_jobs)}) {child_name} [{formula.molecule(child._molecule)}]', ['split']) + log.flow(f'Charge: {child.settings.input.ams.System.charge}', ['straight', 'straight']) + log.flow(f'Spin-Polarization: {child.settings.input.adf.SpinPolarization}', ['straight', 'straight']) + log.flow(f'Work dir: {child.workdir}', ['straight', 'straight']) if child.can_skip(): log.flow(log.Emojis.warning + ' Already ran, skipping', ['straight', 'end']) log.flow() - continue + else: + log.flow(log.Emojis.good + ' Submitting', ['straight', 'end']) + child.run() + self.dependency(child) + log.flow(f'SlurmID: {child.slurm_job_id}', ['straight', 'skip', 'end']) + log.flow() - log.flow(log.Emojis.good + ' Submitting', ['straight', 'end']) - # recast the plams.Settings object into a Result object as that is what run expects - child.settings = results.Result(child_setts[childname]) - child.run() - self.dependency(child) + if self.decompose_elstat: + child_STOFIT = ADFJob(child) + child_STOFIT.name += '_STOFIT' + elstat_jobs[child_STOFIT.name] = child_STOFIT + child_STOFIT.settings.input.adf.STOFIT = '' + child_STOFIT.settings.input.adf.PRINT = 'SFOSiteEnergies Elstat' + child_STOFIT.settings.input.adf.pop("NumericalQuality") + child_STOFIT.settings.input.adf.BeckeGrid.Quality = "Excellent" + + log.flow(f'Fragment ({i}/{len(self.child_jobs)}) {child_name} [{formula.molecule(child._molecule)}] with STOFIT', ['split']) + log.flow(f'Charge: {child_STOFIT.settings.input.ams.System.charge}', ['straight', 'straight']) + log.flow(f'Spin-Polarization: {child_STOFIT.settings.input.adf.SpinPolarization}', ['straight', 'straight']) + log.flow(f'Work dir: {child_STOFIT.workdir}', ['straight', 'straight']) + + if child_STOFIT.can_skip(): + log.flow(log.Emojis.warning + ' Already ran, skipping', ['straight', 'end']) + log.flow() + else: + log.flow(log.Emojis.good + ' Submitting', ['straight', 'end']) + child_STOFIT.run() + self.dependency(child_STOFIT) + log.flow(f'SlurmID: {child_STOFIT.slurm_job_id}', ['straight', 'skip', 'end']) + log.flow() + + child_NoElectrons = ADFJob(child) + child_NoElectrons.name += '_NoElectrons' + elstat_jobs[child_NoElectrons.name] = child_NoElectrons + child_NoElectrons.settings.input.adf.STOFIT = '' + child_NoElectrons.settings.input.adf.PRINT = 'SFOSiteEnergies Elstat' + child_NoElectrons.charge(molecule.number_of_electrons(child_NoElectrons._molecule)) + child_NoElectrons.spin_polarization(0) + child_NoElectrons.settings.input.adf.pop("NumericalQuality") + child_NoElectrons.settings.input.adf.BeckeGrid.Quality = "Excellent" + + log.flow(f'Fragment ({i}/{len(self.child_jobs)}) {child_name} [{formula.molecule(child._molecule)}] without Electrons', ['split']) + log.flow(f'Charge: {child_NoElectrons.settings.input.ams.System.charge}', ['straight', 'straight']) + log.flow(f'Spin-Polarization: {child_NoElectrons.settings.input.adf.SpinPolarization}', ['straight', 'straight']) + log.flow(f'Work dir: {child_NoElectrons.workdir}', ['straight', 'straight']) + + if child_NoElectrons.can_skip(): + log.flow(log.Emojis.warning + ' Already ran, skipping', ['straight', 'end']) + log.flow() + else: + log.flow(log.Emojis.good + ' Submitting', ['straight', 'end']) + child_NoElectrons.run() + self.dependency(child_NoElectrons) + log.flow(f'SlurmID: {child_NoElectrons.slurm_job_id}', ['straight', 'skip', 'end']) + log.flow() - log.flow(f'SlurmID: {child.slurm_job_id}', ['straight', 'skip', 'straight']) - log.flow(f'Work dir: {child.workdir}', ['straight', 'skip', 'end']) - log.flow() # in the parent job the atoms should have the region and adf.f defined as options atom_lines = [] # for each atom we check which child it came from for atom in self._molecule: - for childname, child in self.childjobs.items(): + for child_name, child in self.child_jobs.items(): for childatom in child._molecule: # we check by looking at the symbol and coordinates of the atom if (atom.symbol, atom.x, atom.y, atom.z) == (childatom.symbol, childatom.x, childatom.y, childatom.z): # now write the symbol and coords as a string with the correct suffix - atom_lines.append(f'\t\t{atom.symbol} {atom.x} {atom.y} {atom.z} region={childname} adf.f={childname}') + atom_lines.append(f'\t\t{atom.symbol} {atom.x} {atom.y} {atom.z} region={child_name} adf.f={child_name}') + old_name = self.name # write the atoms block as a string with new line characters self.settings.input.ams.system.atoms = ('\n' + '\n'.join(atom_lines) + '\n\tEnd').expandtabs(4) # set the _molecule to None, otherwise it will overwrite the atoms block self._molecule = None # run this job - self.rundir = j(self.rundir, self.name) + self.rundir = j(self.rundir, old_name) self.name = 'complex' log.flow(log.Emojis.good + ' Submitting parent job', ['split']) super().run() @@ -483,6 +535,7 @@ def run(self): log.flow() # also do the calculation with SCF cycles set to 1 + old_iters = self.settings.input.adf.SCF.Iterations or 300 self.SCF(iterations=0) # we must repopulate the sbatch settings for the new run [self._sbatch.pop(key, None) for key in ['D', 'chdir', 'J', 'job_name', 'o', 'output']] @@ -492,6 +545,74 @@ def run(self): super().run() log.flow(f'SlurmID: {self.slurm_job_id}', ['straight', 'end']) log.flow() + + # also do the calculation with no electrons on the fragments if the user requested a elstat decomposition + if self.decompose_elstat: + frag_names = self.child_jobs.keys() + # reset the SCF iterations + self.SCF(iterations=old_iters) + self.settings.input.adf.pop("NumericalQuality") + self.settings.input.adf.BeckeGrid.Quality = "Excellent" + [self._sbatch.pop(key, None) for key in ['D', 'chdir', 'J', 'job_name', 'o', 'output']] + self.name = 'complex_STOFIT' + + # set the region and fragment files correctly + atom_lines_ = [] + for line in atom_lines: + for frag in frag_names: + line = line.replace(frag, f'{frag}_STOFIT') + atom_lines_.append(line) + + self.settings.input.ams.system.atoms = ('\n' + '\n'.join(atom_lines_) + '\n\tEnd').expandtabs(4) + # set the fragment references correctly + self.settings.input.adf.pop('fragments') + for frag_name in frag_names: + self.settings.input.adf.fragments[frag_name + '_STOFIT'] = j(elstat_jobs['frag_' + frag_name + '_STOFIT'].workdir, 'adf.rkf') + self.settings.input.adf.STOFIT = '' + self.settings.input.adf.PRINT = 'SFOSiteEnergies Elstat' + + log.flow(log.Emojis.good + ' Submitting complex with STOFIT', ['split']) + super().run() + log.flow(f'SlurmID: {self.slurm_job_id}', ['straight', 'end']) + log.flow() + + for frag in frag_names: + # other_frags stores fragment names for fragments that keep their electrons + other_frags = [frag_ for frag_ in frag_names if frag_ != frag] + + # we must repopulate the sbatch settings for the new run + [self._sbatch.pop(key, None) for key in ['D', 'chdir', 'J', 'job_name', 'o', 'output']] + self.name = f'complex_{frag}_NoElectrons' + + # set the region and fragment files correctly + atom_lines_ = [] + for line in atom_lines: + line = line.replace(frag, f'{frag}_NoElectrons') + for other_frag in other_frags: + line = line.replace(other_frag, f'{other_frag}_STOFIT') + atom_lines_.append(line) + self.settings.input.ams.system.atoms = ('\n' + '\n'.join(atom_lines_) + '\n\tEnd').expandtabs(4) + + # set the fragment references correctly + self.settings.input.adf.pop('fragments') + self.settings.input.adf.fragments[frag + '_NoElectrons'] = j(elstat_jobs['frag_' + frag + '_NoElectrons'].workdir, 'adf.rkf') + for other_frag in other_frags: + self.settings.input.adf.fragments[other_frag + '_STOFIT'] = j(elstat_jobs['frag_' + other_frag + '_STOFIT'].workdir, 'adf.rkf') + + other_charge = sum([elstat_jobs['frag_' + other_frag + '_STOFIT'].settings.input.ams.System.charge for other_frag in other_frags]) + total_charge = other_charge + (elstat_jobs['frag_' + frag + '_NoElectrons'].settings.input.ams.System.charge) + self.charge(total_charge) + + other_spin_polarization = sum([elstat_jobs['frag_' + other_frag + '_STOFIT'].settings.input.adf.SpinPolarization for other_frag in other_frags]) + total_spin_polarization = other_spin_polarization + (elstat_jobs['frag_' + frag + '_NoElectrons'].settings.input.adf.SpinPolarization) + self.spin_polarization(total_spin_polarization) + + log.flow(log.Emojis.good + f' Submitting complex with 0 electrons in fragment {frag}', ['split']) + + super().run() + log.flow(f'SlurmID: {self.slurm_job_id}', ['straight', 'end']) + log.flow() + log.flow(log.Emojis.finish + ' Done, bye!', ['startinv']) @@ -588,7 +709,7 @@ def output_cub_paths(self): for sfo in self._sfos: spin_part = '' if sfo.spin == 'AB' else f'_{sfo.spin}' paths.append(f'{cuboutput}%SFO_{sfo.symmetry}{spin_part}%{sfo.index}.cub') - print(paths) + # print(paths) return paths def can_skip(self): diff --git a/src/tcutility/molecule.py b/src/tcutility/molecule.py index 6fbb19cd..67be5479 100644 --- a/src/tcutility/molecule.py +++ b/src/tcutility/molecule.py @@ -5,6 +5,14 @@ from tcutility import ensure_list from tcutility.results import result +from tcutility.data import atom + + +def number_of_electrons(mol: plams.Molecule) -> int: + nel = 0 + for at in mol: + nel += atom.atom_number(at.symbol) + return nel def parse_str(s: str): @@ -88,9 +96,9 @@ def parse_flags(args): for line in atom_lines: # parse every atom first symbol, x, y, z, *args = line.split() - atom = plams.Atom(symbol=symbol, coords=(float(x), float(y), float(z))) - atom.flags = parse_flags(args) - mol.add_atom(atom) + at = plams.Atom(symbol=symbol, coords=(float(x), float(y), float(z))) + at.flags = parse_flags(args) + mol.add_atom(at) # after the atoms we parse the flags for the molecule flag_lines = lines[natoms + 2 :] @@ -189,12 +197,12 @@ def guess_fragments(mol: plams.Molecule) -> Dict[str, plams.Molecule]: return result.Result(fragment_mols) # second method, check if the atoms have a frag= flag defined - fragment_names = set(atom.flags.get("frag") for atom in mol) + fragment_names = set(at.flags.get("frag") for at in mol) if len(fragment_names) > 0: fragment_mols = {name: plams.Molecule() for name in fragment_names} - for atom in mol: + for at in mol: # get the fragment the atom belongs to and add it to the list - fragment_mols[atom.flags.get("frag")].add_atom(atom) + fragment_mols[at.flags.get("frag")].add_atom(at) for frag in fragment_names: fragment_mols[frag].flags = {"tags": set()} @@ -223,9 +231,9 @@ def _xyz_format(mol: plams.Molecule, include_n_atoms: bool = True) -> str: """ n_atoms = len(mol.atoms) if include_n_atoms: - return f"{n_atoms}\n" + "\n".join([f"{atom.symbol:6s}{atom.x:16.8f}{atom.y:16.8f}{atom.z:16.8f}" for atom in mol.atoms]) + return f"{n_atoms}\n" + "\n".join([f"{at.symbol:6s}{at.x:16.8f}{at.y:16.8f}{at.z:16.8f}" for at in mol.atoms]) - return "\n".join([f"{atom.symbol:6s}{atom.x:16.8f}{atom.y:16.8f}{atom.z:16.8f}" for atom in mol.atoms]) + return "\n".join([f"{at.symbol:6s}{at.x:16.8f}{at.y:16.8f}{at.z:16.8f}" for at in mol.atoms]) def _amv_format(mol: plams.Molecule, step: int, energy: Union[float, None] = None) -> str: @@ -240,8 +248,8 @@ def _amv_format(mol: plams.Molecule, step: int, energy: Union[float, None] = Non If no energy is provided, the energy is not included in the string representation""" if energy is None: - return f"Geometry {step}\n" + "\n".join([f"{atom.symbol:6s}{atom.x:16.8f}{atom.y:16.8f}{atom.z:16.8f}" for atom in mol.atoms]) - return f"Geometry {step}, Energy: {energy} Ha\n" + "\n".join([f"{atom.symbol:6s}{atom.x:16.8f}{atom.y:16.8f}{atom.z:16.8f}" for atom in mol.atoms]) + return f"Geometry {step}\n" + "\n".join([f"{at.symbol:6s}{at.x:16.8f}{at.y:16.8f}{at.z:16.8f}" for at in mol.atoms]) + return f"Geometry {step}, Energy: {energy} Ha\n" + "\n".join([f"{at.symbol:6s}{at.x:16.8f}{at.y:16.8f}{at.z:16.8f}" for at in mol.atoms]) def write_mol_to_xyz_file(mols: Union[List[plams.Molecule], plams.Molecule], filename: Union[str, pl.Path], include_n_atoms: bool = False) -> None: @@ -276,9 +284,9 @@ def save(mol: plams.Molecule, path: str, comment: str = None): comment = comment or mol.comment if hasattr(mol, "comment") else "" with open(path, "w+") as f: f.write(f"{len(mol.atoms)}\n{comment}\n") - for atom in mol.atoms: + for at in mol.atoms: flags_str = "" - flags = atom.flags if hasattr(atom, "flags") else {} + flags = at.flags if hasattr(at, "flags") else {} for key, value in flags.items(): if key == "tags": if len(value) > 0: @@ -288,7 +296,7 @@ def save(mol: plams.Molecule, path: str, comment: str = None): else: flags_str += f"{key}={value} " - f.write(f"{atom.symbol}\t{atom.coords[0]: .10f}\t{atom.coords[1]: .10f}\t{atom.coords[2]: .10f}\t{flags_str}\n") + f.write(f"{at.symbol}\t{at.coords[0]: .10f}\t{at.coords[1]: .10f}\t{at.coords[2]: .10f}\t{flags_str}\n") f.write("\n") diff --git a/src/tcutility/results/adf.py b/src/tcutility/results/adf.py index 25ed6af1..aa3e9334 100644 --- a/src/tcutility/results/adf.py +++ b/src/tcutility/results/adf.py @@ -163,7 +163,10 @@ def get_properties(info: Result) -> Result: :Result object containing properties from the ADF calculation: - **energy.bond (float)** – bonding energy (|kcal/mol|). - - **energy.elstat (float)** – total electrostatic potential (|kcal/mol|). + - **energy.elstat.total (float)** – total electrostatic potential (|kcal/mol|). + - **energy.elstat.Vee (float)** – electron-electron repulsive term of the electrostatic potential (|kcal/mol|). + - **energy.elstat.Ven (float)** – electron-nucleus attractive term of the electrostatic potential (|kcal/mol|). + - **energy.elstat.Vnn (float)** – nucleus-nucleus repulsive term of the electrostatic potential (|kcal/mol|). - **energy.orbint.total (float)** – total orbital interaction energy containing contributions from each symmetry label and correction energy(|kcal/mol|). - **energy.orbint.{symmetry label} (float)** – orbital interaction energy from a specific symmetry label (|kcal/mol|). - **energy.orbint.correction (float)** - orbital interaction correction energy, the difference between the total and the sum of the symmetrized interaction energies (|kcal/mol|) @@ -198,7 +201,29 @@ def get_properties(info: Result) -> Result: # read energies (given in Ha in rkf files) ret.energy.bond = reader_adf.read("Energy", "Bond Energy") * constants.HA2KCALMOL - ret.energy.elstat = reader_adf.read("Energy", "elstat") * constants.HA2KCALMOL + + # total electrostatic potential + ret.energy.elstat.total = reader_adf.read("Energy", "elstat") * constants.HA2KCALMOL + + # we can further decompose elstat if it was enabled + if info.files.out: + with open(info.files.out) as output: + lines = output.readlines() + + skip_next = -1 + for line in lines: + if "Electrostatic Interaction Energies" in line: + skip_next = 4 + continue + if skip_next == 0: + f1, f2, Vee, Ven, Vnn, total = line.strip().split() + ret.energy.elstat.Vee = float(Vee) * constants.HA2KCALMOL + ret.energy.elstat.Ven = float(Ven) * constants.HA2KCALMOL + ret.energy.elstat.Vnn = float(Vnn) * constants.HA2KCALMOL + skip_next -= 1 + + + # print(info.files) # read the total orbital interaction energy diff --git a/src/tcutility/results/ams.py b/src/tcutility/results/ams.py index 6b7e7786..102d2be6 100644 --- a/src/tcutility/results/ams.py +++ b/src/tcutility/results/ams.py @@ -71,6 +71,9 @@ def get_calc_files(calc_dir: str) -> dict: f = ".".join(parts[1:]).replace("logfile", "log") ret[f] = os.path.abspath(file) + if file.endswith(".out") and "CreateAtoms" not in file: + ret["out"] = os.path.abspath(file) + return ret @@ -274,10 +277,7 @@ def _get_fragment_indices_from_input_order(results_type) -> arrays.Array1D: """Function to get the fragment indices from the input order. This is needed because the fragment indices are stored in internal order in the rkf file.""" frag_indices = np.array(results_type.read("Geometry", "fragment and atomtype index")).reshape(2, -1) # 1st row: Fragment index; 2nd row atomtype index atom_order = np.array(results_type.read("Geometry", "atom order index")).reshape(2, -1) # 1st row: input order; 2nd row: internal order - combined = np.concatenate((frag_indices, atom_order), axis=0) - sorted_by_input_order = np.array(sorted(combined.T, key=lambda x: x[2])) # sort the fragment indices by input order - frag_order = sorted_by_input_order[:, 0] - return frag_order + return frag_indices[0][atom_order[0] - 1] # ------------------------ Molecule -------------------------- # diff --git a/test/test_vdd_analysis.py b/test/test_vdd_analysis.py index 48f0e190..4caf40b9 100644 --- a/test/test_vdd_analysis.py +++ b/test/test_vdd_analysis.py @@ -65,7 +65,7 @@ def vdd_manager_geo_nosym(): def test_get_fragment_indices_from_input_order_disordered(kfreader_fa_disordered_fragindices_cs): frag_indices = ams._get_fragment_indices_from_input_order(kfreader_fa_disordered_fragindices_cs) - assert np.array_equal(frag_indices, np.array([1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1])) + assert np.array_equal(frag_indices, np.array([1, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2])) def test_get_fragment_indices_from_input_order_ordered(kfreader_fa_ordered_fragindices_cs):