duartegroup / autode Goto Github PK

Describe the bug

I.
When a Hessian matrix is to be calculated, it does not work.
The error message is CouldNotGetProperty: Could not get the Hessian: could not convert string to float: '1\H'.
I tried with other SMILEs string and they still do no work.

II.
The bottom line block (i.e. 1\1\ ... @\n) of log file generated after running Gaussian16 has "3 * number_of_atoms". In this case 6. Due to this, ValueError: Array was not the correct shape to be broadcast into the lower triangle. Need N(N+1)/2 elements for an NxN array is made.
Is the bottom block really hessian information? (I am not versed in QM calculation).

To Reproduce

import autode
h2 = autode.Molecule(smiles='[H][H]')
calc=autode.Calculation(name='h2',molecule=h2,method=autode.wrappers.G16.G16(),keywords=autode.wrappers.G16.G16().keywords.opt)
calc.run()
calc.get_energy() #Works!
calc.get_gradients() #Works!
calc.get_hessian() #Breaks with msg: CouldNotGetProperty: Could not get the Hessian: could not convert string to float: '1\H'

Expected behavior

Hessian matrix should be calculated.

Environment

• OS: Centos8
• Version: 4.18.0-147.5.1.el8_1.x86_64

Suggestion
I.
In autode/wrappers/G09.py file,

CHANGE from
hess_str = "".join(hess_lines[::-1]).split(r'\')[-3]
hess_values = [float(val) for val in hess_str.split(',')]
TO
hess_str = list(map(lambda inp: inp.split("\")[0], hess_str.split(',')[2:])) #CouldNotGetProperty error will be fixed as a result of removing \H, \C character strings, in addition to removing '1' and '1' characters.
hess_values = [float(val) for val in hess_str.split(',')] #This will save values correctly into a float value list.

II.
Is the bottom line block (i.e. 1\1\ ... @\n) really Hessian? In this case 6 elements (3 * num_atoms) exist. However, this information seems to be exactly the same as optimized Cartesian coordinate (not Hessian!), accessible by calc.get_final_atoms(). If so, this should be corrected and correct method to calculate (or read) Hessian matrix should be updated.

Additional context

Hello again!

I tried using PBEh-3c as optimisation method and that seems to work great so far. They suggest, however, to use a scaling factor of 0.95 for computing thermodynamics. I added the following bit to the hess-keywords which works fine.

%freq
scalfreq 0.95
end

ORCA 5.0.3 produces the following after autodE called for the calculation.

-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  0.950000000  (already applied!)

   0:         0.00 cm**-1
   1:         0.00 cm**-1
   2:         0.00 cm**-1
   3:         0.00 cm**-1
   4:         0.00 cm**-1
   5:         0.00 cm**-1
   6:      1620.46 cm**-1
   7:      3708.47 cm**-1
   8:      3815.22 cm**-1

Checking the frequencies in autodE, however, gives:

if molecule.imaginary_frequencies:
    frequencies = molecule.imaginary_frequencies + molecule.vib_frequencies
else:
    frequencies = molecule.vib_frequencies
print(*[f'{(freq).to("cm-1"):.1f}' for freq in frequencies],sep=', ')

Frequencies in cm-1
1705.8, 3903.8, 4016.2

This is weird, because I cannot find the numbers at all in ORCA's output files. Does autodE re-compute the frequencies straight from the found normal modes? Both in the .hess and the .out, ORCA only stores the scaled frequencies.

Do I assume correctly that also the thermodynamics did not get the scaling factor applied? My test computing water several times seems to support that. The deviations appear to be in our margin of error (full cycles with conformer search and opt were run each time).

For H2O:

No scaling
E_opt, G_cont, H_cont, E_sp
-76.267818102248,0.0073088101511278696,0.0257081115702916,-76.385195142966

Strong scaling (0.5)
E_opt, G_cont, H_cont, E_sp
-76.26781810866,0.007308818432288296,0.025708124156216056,-76.385195140162
-76.26781810866,0.0073088378537382,0.02570814342764412,-76.385195140162
-76.267818106119,0.007308850889597047,0.025708153656751487,-76.385195141268
-76.267818102248,0.007308814161775101,0.02570811555037913,-76.385195142966
-76.26781810866,0.007308829535297789,0.025708135175725795,-76.385195140162

Is there a good way to transfer the scaling factor then (and propagate it to the stored H_cont and G_cont)? I could surely overwrite the frequencies with the scaled ones but how would I go best about re-computing H_cont and G_cont without triggering new frequency calculations?

Or better: can you add a snippet to ade.Config to define a scaling factor that is given before the thermodynamics are computed?

Describe the bug

When I input:
Nu = ade.Reactant(smiles='[S-]C', charge = -1)
It gives me error:
UFFTYPER: Unrecognized charge state for atom: 0

To Reproduce

Nu = ade.Reactant(smiles='[S-]C', charge = -1)

Expected behavior

I expect methyl thiolate anion can be recognized

Environment

• Operating System: Red Hat Linux
• Python version: 3.8.2
• autodE version: 1.2.2

Additional context

Hi there, is it possible for the module to generate a z-matrix from the xyz coordinates that retains the atom number labelling and follows a carbon skeleton. I would like to use the module to build a starting transition state geometry by adding on a reaction centre and relating the position to a given atom number in the compound structure and then feed that into a transition state search. Is it then possible to convert it back to xyz coordinates again conserving atom number labelling?

Many thanks

We have a quick question:

Can one use autoDE with a PES from a ML method, such as gap, or ani ?

Do you have an example ?

in autoDE/wrappers/methods.py there are some classes like this one:

class ExternalMethodEGH(ExternalMethod, ABC):

which I guess is what we need, but better if you had an example.

Thanks !

Describe the bug
When running an XTB geometry optimisation, a CouldNotGetProperty(energy) error is raised, despite the total energy being available in the XTB output file. Looking at the source-code in XTB.py (191-199) the line containing the total energy in the output file is expected in autodE to be in capitals (i.e. TOTAL ENERGY), but in the .out file is lowercase 'total energy'. Changing this line locally on the source code fixes the big.

To Reproduce

mol = autode.Molecule(name='test', smiles='FC1=CC=C(NC(N2CCCC2)=O)C=C1C3=NC4=NC=C(N5CCOCC5)N4C=C3')
mol.optimise(method=ade.methods.XTB())

This raises a CouldNotGetProperty error

Expected behavior
The energy to be located, and the optimisation to complete normally

Environment

• Operating System: MacOS Monterey 12.2.1
• Python version: 3.9.6
• autodE version: 1.2.1
• XTB version: 6.4.1

Additional context

Hello again,

I ran into a reaction class that should be easy enough to implement but does not seem to be implemented at the moment:

C/C=C/C>>C/C=C\C

E/Z-isomerisation by double bond rotation. I am sure with the current framework that I can work around it, but the reaction optimiser should perhaps learn that?

All the best,
Florian

Describe the bug

Although having a structure with >9 rings is uncommon it is present in some of the autodE example metal complexes shown in the paper – the new SMILES parser should be able to parse them.

To Reproduce

import autode as ade
mol = ade.Molecule(smiles='[N+]12=CC=CC3=C1C(C(C=C3)=CC=C4)=[N+]4[Cr]256([N+]7=CC=CC8=C7C9=[N+]6C=CC=C9C=C8)[N+]%10=CC=CC%11=C%10C%12=[N+]5C=CC=C%12C=C%11')
...
  File "/Users/tom/miniconda3/envs/test/lib/python3.9/site-packages/autode/smiles/parser.py", line 76, in smiles
    self._check_smiles()
  File "/Users/tom/miniconda3/envs/test/lib/python3.9/site-packages/autode/smiles/parser.py", line 63, in _check_smiles
    raise InvalidSmilesString(f'{self._string} had invalid characters:'
autode.exceptions.InvalidSmilesString: ... had invalid characters:['%']
>>> 

Expected behavior

SMILES should be able to be parsed and the structure built

Environment

• OS: Mac OS
• Version: 1.0.3_6

Additional context

dear autodE developers,
I attach a patch for automatically adding ecps in nwchem input under the def2 basis set.
Thank you for your work
fulvio

autode_nwchem_diff.txt

Hi there,

I just have a tiny question. I will want to automatise (e.g. via RDKit) estimating a good number of cores to request for autodE calcs. I saw that you preferably send many 1-core jobs (smart). This + different amounts of active degrees of freedom means that the number of electrons will not be a good sole indicator for this. I do not want to go overboard with the number of cores and just always request the maximum available, as there will be cases where we end up with small molecules and only a few conformers. Here we would then probably still end up putting many more cores than needed on tiny molecules. So some sort of upper limit would be good as best practice to keep the efficiency high, right?

Would you happen to have a good measure for guesstimating a good amount of cores or did you even run some benchmarks? I would guess something like 3^number_of_rotatable_bonds x (numbers_of_electrons)^3 x factor_of_how_rich_my_cluster_is should be decent? I will add an upper limit and perhaps a correction for rotatable bonds that do not actually change anything (e.g. -CH3).

Thank you!

Possible API Improvement

Currently with a Calculation calc.get_energy() returns None if a calculation failed but raises a CalculationError for e.g. calc.get_gradients(); this is not particularly consistent or user friendly.

Having talked it over with @tristan-j-wood having properties (e.g. calc.energy) that are set after calc.run(), and may be None (after throwing a logging error), would be easier and more consistent with the rest of the API. Also having access to the calc.molecule properties directly e.g.

class Molecule:

    def __init__(self, x):

        self.x = x


class Calculation:

    def __getattr__(self, item):
        if item in self.molecule.__dict__:
            return self.molecule.__dict__[item]

        return self.__dict__[item]

    def __init__(self, molecule):

        self.molecule = molecule


if __name__ == '__main__':

    calc = Calculation(molecule=Molecule(x=1))

    print(calc.x)  # -> 1

    calc.x = 2
    print(calc.x)  # -> 2

would be more readable. Also should throw more useful exceptions if calculations are initialised but not run before getting a property.

Describe the bug
WHen performing TS search, calculation stops with that message:
Could not build a molecular graph

Please check attachment for reaction indication

Thanks a lot!!

To Reproduce
Python 3.9.10 | packaged by conda-forge | (main, Jan 30 2022, 18:04:04)
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

import autode as ade
ade.Config.hcode = 'NWChem'
ade.Config.lcode = 'XTB'
ade.Config.n_cores = 2
ade.methods.get_lmethod()
XTB(available = True)
ade.methods.get_hmethod()
NWChem(available = True)
rxn = ade.Reaction('CN1OC(N+=O)CC1>>CN2OC=CC2.O=NO', name='PIR1')
rxn.calculate_reaction_profile()
Traceback (most recent call last):
File "", line 1, in
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/utils.py", line 344, in wrapped_function
return func(*args, **kwargs)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/reactions/reaction.py", line 138, in calculate_reaction_profile
calculate(self)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/utils.py", line 162, in wrapped_function
result = func(*args, **kwargs)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/reactions/reaction.py", line 126, in calculate
reaction.find_lowest_energy_conformers()
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/reactions/reaction.py", line 511, in find_lowest_energy_conformers
mol.find_lowest_energy_conformer(hmethod=h_method)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/utils.py", line 162, in wrapped_function
result = func(*args, **kwargs)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/species/species.py", line 1228, in find_lowest_energy_conformer
self.conformers.prune_diff_graph(self.graph)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/conformers/conformers.py", line 202, in prune_diff_graph
make_graph(conformer)
File "/home/cnieto/.conda/envs/autode_env/lib/python3.9/site-packages/autode/mol_graphs.py", line 83, in make_graph
raise ex.NoAtomsInMolecule('Could not build a molecular graph with no '
autode.exceptions.NoAtomsInMolecule: Could not build a molecular graph with no atoms

Expected behavior

Environment

• Operating System: Ubutnu 21.04 LTS
• Python version: 3.9.10
• autodE version: 1.2.1

Additional context

The document says the Config.max_core sets the memory per core. But using Gaussian16 (didn't test g09 yet), the %mem (total memory use) in the input file is set to the same value as Config.max_core. This might leads to insufficient memory use that causes errors in some clusters or just slows down the job.

Describe the bug
SMILES parser generates the wrong Multiplicity

To Reproduce
import autode as ade

ade.Config.n_cores = 28
ade.Config.hcode = 'g16'

ade.Config.G16.keywords.set_opt_basis_set('6-31+G(d)')
ade.Config.G16.keywords.sp.basis_set = '6-311+G(d,p)'
ade.Config.G16.keywords.set_opt_functional('B3LYP')
ade.Config.G16.keywords.set_functional('B3LYP')

test = ade.Molecule(name='test',smiles='CCC1=CC=C(C=C1)C[Na]')
test.print_xyz_file()

Expected behavior
The expected results are as follows
INFO Initialisation with SMILES successful. Charge=0, Multiplicity=1, Num. Atoms=21

Environment

• OS: Unix
• Version: 1.0.3

Additional context

INFO Generating a Molecule object for test
INFO Parsing CCC1=CC=C(C=C1)C[Na]
INFO Inverting stereochemistry
INFO Parsed SMILES in: 0.42 ms
INFO Setting 21 atom types
INFO Have 1 ring(s)
INFO Queue: [1]
INFO Queue: [2]
INFO Queue: [7, 3]
INFO Queue: [3, 6]
INFO Queuing Dihedral(idxs=[1, 2, 3, 4], φ0=3.14)
INFO Have 1 dihedral(s) to rotate
INFO Resetting sites on atom 4
INFO Rotating 3 sites onto 1 points
INFO Resetting sites on atom 6
INFO Rotating 3 sites onto 1 points
INFO Queue: [6, 4]
INFO Queuing Dihedral(idxs=[5, 6, 7, 2], φ0=0)
INFO Have 1 dihedral(s) to rotate
WARNING Could not apply rotation Dihedral(idxs=[5, 6, 7, 2], φ0=0)
INFO Queue: [4, 5]
INFO Closing ring on: SMILESBond([4, 5], order=2) and adjusting atoms
INFO Adjusting ring dihedrals to close the ring
INFO Closed ring in: 1.66 ms
INFO A ring was poorly closed - adjusting angles
WARNING Closing large rings not implemented
WARNING Failed to close a ring, minimising on all atoms
INFO Double bond - adding constraint
INFO Double bond - adding constraint
INFO Double bond - adding constraint
INFO Resetting sites on atom 4
INFO Rotating 3 sites onto 2 points
INFO Resetting sites on atom 5
INFO Rotating 3 sites onto 2 points
INFO Queue: [5]
INFO Queuing Dihedral(idxs=[3, 4, 5, 6], φ0=0)
INFO Have 1 dihedral(s) to rotate
WARNING Could not apply rotation Dihedral(idxs=[3, 4, 5, 6], φ0=0)
INFO Queue: [8]
INFO Queue: []
INFO Minimising non-bonded repulsion by dihedral rotation
INFO Have 1 dihedrals to rotate
INFO Performed final dihedral rotation in: 9.02 ms
INFO Built 3D in: 63.47 ms
INFO Generating molecular graph with NetworkX
INFO Generating molecular graph with NetworkX
INFO Setting the π bonds in a species
INFO Setting the stereocentres in a species
WARNING Could not find a valid valance for Na. Guessing at 6
INFO Initialisation with SMILES successful. Charge=0, Multiplicity=2, Num. Atoms=21

Describe the bug

SMILES parser generates the wrong connectivity for ring closures following branches

To Reproduce

from autode import Molecule
from autode.smiles.parser import Parser
from autode.smiles.builder import Builder

parser = Parser()
builder = Builder()
parser.parse('CC12[C@@]3(CCC4)C4=C[C@@H](C1C=CO2)S(=O)3=O')
builder.build(parser.atoms, parser.bonds)

mol = Molecule(atoms=builder.atoms)
mol.print_xyz_file(filename='tmp.xyz')

Expected behavior

SO2 should form a 5-membered ring not a 6-membered with the O part of the ring

Environment

• OS: Mac OS
• Version: 1.0.3

Additional context

Describe the bug
DA test run ended with abnormal TS profile.

Any hints of the source for this?

To Reproduce

Python 3.9.10 | packaged by conda-forge | (main, Jan 30 2022, 18:04:04)
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

import autode as ade
ade.Config.hcode = 'NWChem'
ade.methods.get_hmethod()
NWChem(available = True)
ade.methods.get_lmethod()
XTB(available = True)
ade.Config.n_cores = 8
rxn = ade.Reaction('C=CC=C.C=C>>C1=CCCCC1', name='DA')
rxn.calculate_reaction_profile()

Expected behavior

Environment

• Operating System: Linux, Ubuntu 20.04 LTS
• Python version:
• autodE version:

Additional context

Describe the bug
When I go to import autode it produces the following error message
ImportError: /home/maihfly/anaconda3/envs/Cyclodextrins-env-2_0/lib/python3.9/site-packages/rdkit/Chem/../../../../libRDKitmaeparser.so.1: undefined symbol: _ZN5boost9iostreams4zlib8deflatedE
To Reproduce
you import autode as ade and this error appears

Expected behavior
I was expecting autode to import

Environment

• Operating System: lastes ubuntu linux os at the time of writing
• Python version: 3.9
• autodE version: 1.2.0

Additional context

Describe the bug

NCI complex conformers are all generated with the same name, which means optimisations skip over all but the first one. The documentation water trimer works, but for other more complex surfaces this is not expected to.

To Reproduce

import autode as ade
from autode.species import NCIComplex

ade.Config.num_complex_sphere_points = 5
ade.Config.num_complex_random_rotations = 3

water = ade.Molecule(name='water', smiles='O')

trimer = NCIComplex(water, water, water, name='water_trimer')
trimer.find_lowest_energy_conformer(lmethod=ade.methods.XTB())
trimer.print_xyz_file()

The conformers directory only contains a single XTB output file

Expected behavior

Environment

• OS: CentOS
• Version: v. 1.1.0

Hi, I have performed the following test and it failed:

import autode as Ade
mol = ade.Molecule("water.xyz")
python
Python 3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:41:54) [Clang 13.0.1 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import autode as ade
>>> mol = ade.Molecule("example.xyz")
>>> xtb = ade.methods.XTB()
>>> calc = ade.Calculation(name="water", molecule=mol, method=xtb, keywords=xtb.keywords.hess)
>>> mol.calc_thermo(calc=calc)
Traceback (most recent call last):
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/calculation.py", line 413, in get_hessian
    hessian = Hessian(self.method.get_hessian(self),
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/utils.py", line 356, in wrapped_function
    assert hasattr(args[0], 'output')
AssertionError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/utils.py", line 285, in wrapped_function
    return func(*args, **kwargs)
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/species/species.py", line 1126, in calc_thermo
    self._run_hess_calculation(method=method, calc=calc, keywords=keywords)
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/species/species.py", line 797, in _run_hess_calculation
    self.hessian = calc.get_hessian()
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/calculation.py", line 38, in wrapped_function
    return func(*args, **kwargs)
  File "/Users/leticiamadureira/opt/anaconda3/envs/kinetic/lib/python3.10/site-packages/autode/calculation.py", line 421, in get_hessian
    raise ex.CouldNotGetProperty(f'Could not get the Hessian: {err}')
autode.exceptions.CouldNotGetProperty: Could not get the Hessian:
>>>

I use you example ,but it still wrong ,why?

from autode import Product
from autode import Reactant
from autode import Reaction
r1 = Reactant('2.xyz',charge=1,mult=3)
r2 = Reactant('3.xyz')
r3 = Product('4.xyz',charge=1,mult=3)
reaction = Reaction(r1, r3,r2, name='cope')
reaction.locate_transition_state()
print(reaction.tss())
Traceback (most recent call last):
File "", line 1, in
TypeError: 'list' object is not callable
print(reaction.tss)
[<autode.transition_states.transition_state.TransitionState object at 0x1519ea8f2910>]
reaction.tss()
Traceback (most recent call last):
File "", line 1, in
TypeError: 'list' object is not callable
a=reaction.tss
print(a)
[<autode.transition_states.transition_state.TransitionState object at 0x1519ea8f2910>]
print(a())
Traceback (most recent call last):
File "", line 1, in
TypeError: 'list' object is not callable
a=reaction.tss()
Traceback (most recent call last):
File "", line 1, in
TypeError: 'list' object is not callable

Describe the bug

Multiplicities defined in ade.Molecule constructors are not always respected.

To Reproduce

>>> import autode as ade
>>> mol = ade.Molecule(smiles='[Sc]C', mult=3)
>>> mol.mult
1

Expected behavior

mol.mult should be 3

Environment

• OS: Linux
• Version: 1.0.5

Hello! Please, excuse my ignorance, but I am trying to figure out how to change number of geom_opt_max_cycles for any type of calculation with QChem, but all I was able to to do is:

import autode as ade
from autode.input_output import xyz_file_to_atoms
from autode.wrappers.keywords import MaxOptCycles

xtb = ade.methods.XTB()
qchem = ade.methods.QChem()
qchem.keywords.low_opt = ['method pbe0', 'basis def2-SVP', 'jobtype opt', MaxOptCycles(15), 'dft_d D3_BJ']

this is for
for x_water_trimer.find_lowest_energy_conformer(lmethod=xtb, hmethod=qchem, allow_connectivity_changes=False)

But in here I overwrite all keywords. If I specify only MaxOptCycles(15), other keywords are omitted and calculation will not run at all. What is the most convenient and correct way for this?
I am trying to di it now for non covalent interaction complexes and for reaction calculation. Any advice will be very appreciated!
Thank you!

Describe the bug
get_atomic_charges() returns a list of 0s for ORCA calculations, when atomic charges are present in the .out file

To Reproduce

mol = autode.Molecule(smiles='c1ncccc1', name='pyridine')
calc = autode.Calculation(name='test_calculation', method=autode.methods.ORCA(), keywords=autode.Config.ORCA.keywords.opt, molecule=mol)
calc.run()
print(calc.get_atomic_charges())

This will print a list of 0s. When the same is repeated with XTB a list of sensible atomic charges is returned.

Expected behavior
A list of non 0 atomic charges should be returned

Environment

• Operating System: MacOS Monterey 12.2.1
• Python version: 3.10.4
• autodE version: 1.2.3
• ORCA version: 4.2.1

Travis CI builds have ceased on travis-ci.org so this is a good opportunity to migrate to GitHub Actions. Add:

• Build/test: https://github.com/t-young31/gap-train/blob/master/.github/workflows/pytest.yml
• CodeCov: https://github.com/codecov/codecov-action
• Deploy: https://github.com/marketplace/actions/deploy-to-github-pages

Hello,

I managed to install autodE using conda (conda install autode --channel conda-forge) and when I ran the tests, all 397 passed. However, when I'm trying the reactions given under the example I'm getting the following errors:

diels_alder - File "/autode/conformers/conformers.py", line 173, in _parallel_calc
n_cores_pp = max(Config.n_cores // len(self), 1)
ZeroDivisionError: integer division or modulo by zero

sn2 - File "/autode/wrappers/ORCA.py", line 499, in get_hessian
raise CouldNotGetProperty('Could not find Hessian file')
autode.exceptions.CouldNotGetProperty: Could not find Hessian file

elimination - File "/autode/wrappers/XTB.py", line 284, in get_final_atoms
raise AtomsNotFound
autode.exceptions.AtomsNotFound

What am I missing here?

Describe the bug

dear friends,could you give me the example of find_tss, I still can not write it correctly. And another method also does not give a correct result
To Reproduce

import autode as ade
ade.Config.n_cores =1
rxn = ade.Reaction('C=CC=C.C=C>>C1=CCCCC1', name='DA')
rxn.locate_transition_state()
print(rxn.tss)
[<autode.transition_states.transition_state.TransitionState object at 0x1542af0c44f0>]

（1） locate_transition_state only locates a single transtion state for each possible bond rearrangment and does not attempt to search the conformational space. how can I find more transtion state for each possible bond rearrangment .
(2) For a reaction step, if there are multiple transition states, how to find them?
（3）can find_tss achieve my goal?
(4)hope you can give me a completel example，thanks very much
Expected behavior

Environment

• OS:
• Version:

Additional context

Hi there,

I assume you guys will use some form of queue manager for running autode calculations. I am struggling to get autode running in SLURM with ORCA as ORCA has this weird behaviour of wanting to run mpiexec itself. Giving SLURM sufficient cpus-per-task crashes the calculation. Giving it sufficient tasks (which is how running ORCA itself works) just runs autode several times. Does anyone have some pointers on this one or even an example SLURM script?

Thank you!
Florian

I recently installed autodE and was trying to reproduce the butadiene-ethene DA reaction, using python 3.7.4 and orca 5.0.2 (although the same error was reproduced using orca 4.1.0). However, the calculation fails during the transition state search phase with the following error:

Process Process-143:
Traceback (most recent call last):
File "/cluster/apps/nss/python/3.7.4/x86_64/lib64/python3.7/multiprocessing/process.py", line 297, in _bootstrap
self.run()
File "/cluster/apps/nss/python/3.7.4/x86_64/lib64/python3.7/multiprocessing/process.py", line 99, in run
self._target(*self._args, **self._kwargs)
File "/cluster/home/rueppf/.local/lib/python3.7/site-packages/autode/utils.py", line 380, in handler
queue.put(func(*args, **kwargs))
File "/cluster/home/rueppf/.local/lib/python3.7/site-packages/autode/mol_graphs.py", line 484, in is_isomorphic
graph1, graph2 = get_graphs_ignoring_active_edges(graph1, graph2)
File "/cluster/home/rueppf/.local/lib/python3.7/site-packages/autode/mol_graphs.py", line 451, in get_graphs_ignoring_active_edges
for edge in ga.edges:
File "/cluster/apps/nss/python/3.7.4/x86_64/lib64/python3.7/site-packages/networkx/classes/reportviews.py", line 1028, in __iter__
for nbr in nbrs:
RuntimeError: dictionary changed size during iteration

The code to obtain this result is simply:

import autode as ade
ade.Config.n_cores = 8
ade.Config.max_core = 4000
rxn = ade.Reaction('C=CC=C.C=C>>C1=CCCCC1', name='DA')
rxn.calculate_reaction_profile()

Environment

• Python version: 3.7.4
• autodE version: 1.2.2
• ORCA version: 5.0.2 (also with 4.1.0)

Hi!

I'm already stressed Duarte group gonna hate me and @t-young31 in particular, lol, but:

I'm working on acetylene behavior in water for CuAAC reaction and conformer search works great as it allows to get more good structures rather than me guessing where to put water molecules. But I am now moving to the next step when I want to use Cu(I) coordinated to triple bond and find lowest complex with several water molecules around. I use an xyz input for alkyne-Cu(I) complex. Unfortunately, autodE gives me the following error in the INFO file:

2022-04-08 21:40:11 autode.log.log: INFO     Getting atoms from PhC2H_CuCl.xyz
2022-04-08 21:40:11 autode.log.log: INFO     Generating a Molecule object for molecule
2022-04-08 21:40:11 autode.log.log: WARNING  Not requested any solvent - returning None
2022-04-08 21:40:11 autode.log.log: INFO     Generating molecular graph with NetworkX
2022-04-08 21:40:11 autode.log.log: INFO     Setting the π bonds in a species
2022-04-08 21:40:11 autode.log.log: WARNING  Cu not found in π valency dictionary - assuming not a π-atom
2022-04-08 21:40:11 autode.log.log: WARNING  Cu not found in π valency dictionary - assuming not a π-atom
2022-04-08 21:40:11 autode.log.log: WARNING  Cu not found in π valency dictionary - assuming not a π-atom
2022-04-08 21:40:11 autode.log.log: INFO     Setting the stereocentres in a species
2022-04-08 21:40:11 autode.log.log: WARNING  Could not find a valid valance for Atom(Cu, 0.7681, -0.3143, -2.5625). Guessing at 6

The ade input I use is basically similar to water trimer example and it worked great for a number of alkynes I used so far.

Am I missing something? Any suggestions would be awesome, thank you!

Description

Nudged elastic band (NEB) interpolation is performed linearly, which can result in non-ideal geometries with small/large bond lengths which hinder convergence.

To Reproduce

For example, rotating a methane molecule and interpolating

import autode as ade

neb = ade.neb.NEB(initial_species=ade.Molecule('methane.xyz'),
                  final_species=ade.Molecule('methane_rot.xyz'),
                  num=10)
neb.interpolate_geometries()

combined = ade.Molecule()
for image in neb.images:
    combined.atoms += image.species.atoms

combined.print_xyz_file(filename='interpolated.xyz')

generates these intermediate geometries

Suggested Fix

Implement the IDPP method from this paper. Or, at the very least, throw an exception if a NEB is calculated from an initial set of geometries that aren't reasonable.

Environment

• OS: MacOS
• autodE version: 1.1.3

Files to reproduce:
data.zip

Describe the bug

When running the reaction below, it crashes with the following:

NoConformers Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 rxn.calculate_reaction_profile()

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/utils.py:345, in requires_hl_level_methods..wrapped_function(*args, **kwargs)
340 except MethodUnavailable:
341 raise MethodUnavailable(f'Function {func.name} requires both'
342 f' a high and low-level method but '
343 f'{suffix}')
--> 345 return func(*args, **kwargs)

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/reactions/reaction.py:138, in Reaction.calculate_reaction_profile(self, units, with_complexes, free_energy, enthalpy)
135 reaction.print_output()
136 return None
--> 138 calculate(self)
140 if not with_complexes:
141 plot_reaction_profile([self], units=units, name=self.name,
142 free_energy=free_energy, enthalpy=enthalpy)

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/utils.py:163, in work_in..func_decorator..wrapped_function(*args, **kwargs)
160 os.mkdir(dir_path)
162 os.chdir(dir_path)
--> 163 result = func(*args, **kwargs)
164 os.chdir(here)
166 if len(os.listdir(dir_path)) == 0:

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/reactions/reaction.py:126, in Reaction.calculate_reaction_profile..calculate(reaction)
124 @work_in(self.name)
125 def calculate(reaction):
--> 126 reaction.find_lowest_energy_conformers()
127 reaction.optimise_reacs_prods()
128 reaction.locate_transition_state()

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/reactions/reaction.py:511, in Reaction.find_lowest_energy_conformers(self)
508 h_method = get_hmethod() if Config.hmethod_conformers else None
509 for mol in self.reacs + self.prods:
510 # .find_lowest_energy_conformer works in conformers/
--> 511 mol.find_lowest_energy_conformer(hmethod=h_method)
513 return None

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/utils.py:163, in work_in..func_decorator..wrapped_function(*args, **kwargs)
160 os.mkdir(dir_path)
162 os.chdir(dir_path)
--> 163 result = func(*args, **kwargs)
164 os.chdir(here)
166 if len(os.listdir(dir_path)) == 0:

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/species/species.py:1230, in Species.find_lowest_energy_conformer(self, lmethod, hmethod, allow_connectivity_changes)
1227 if not allow_connectivity_changes:
1228 self.conformers.prune_diff_graph(self.graph)
-> 1230 self._set_lowest_energy_conformer()
1231 logger.info(f'Lowest energy conformer found. E = {self.energy}')
1232 return None

File ~/.conda/envs/autode_env/lib/python3.10/site-packages/autode/utils.py:317, in requires_conformers..wrapped_function(*args, **kwargs)
314 assert hasattr(args[0], 'n_conformers')
316 if args[0].n_conformers == 0:
--> 317 raise NoConformers
319 return func(*args, **kwargs)

NoConformers:

To Reproduce

Expected behavior

Any hints about what is going on? NWChem conformer opt finished ok.

Environment

• Operating System: Ubutnu 21.04 LTS
• Python version: Python 3.10.4
• autodE version: latest

Additional context

Describe the bug
Once installation and proper H2 exercise was succesfully completed, next example (DA) was performed.
Failed to reproduce it at the light of "autode.exceptions.NoConformers"

To Reproduce
Python 3.9.10 | packaged by conda-forge | (main, Jan 30 2022, 18:04:04)
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

import autode as ade
ade.methods.get_lmethod()
MOPAC(available = True)
ade.Config.hcode = 'NWChem'
ade.methods.get_hmethod()
NWChem(available = True)
ade.Config.n_cores = 8
rxn.calculate_reaction_profile()
Traceback (most recent call last):
File "", line 1, in
File "/home/user/autodE/autode/reactions/reaction.py", line 566, in calculate_reaction_profile
calculate(self)
File "/home/user/autodE/autode/utils.py", line 116, in wrapped_function
result = func(*args, **kwargs)
File "/home/user/autodE/autode/reactions/reaction.py", line 554, in calculate
reaction.find_lowest_energy_conformers()
File "/home/user/autodE/autode/reactions/reaction.py", line 340, in find_lowest_energy_conformers
mol.find_lowest_energy_conformer(hmethod=h_method)
File "/home/user/autodE/autode/utils.py", line 116, in wrapped_function
result = func(*args, **kwargs)
File "/home/user/autodE/autode/species/species.py", line 988, in find_lowest_energy_conformer
self._set_lowest_energy_conformer()
File "/home/user/autodE/autode/utils.py", line 259, in wrapped_function
raise NoConformers
autode.exceptions.NoConformers

Expected behavior
Successful completion of DA example

Environment

• OS: Ubuntu
• Version: 20.04.3 LTS

Additional context

Describe the bug
After installing AutoDe under a conda environment, its unable to find NWChem installation.
It is already in the PATH of my ~/.bashrc file and additionally I have to introduce it at config.py file

path ='/apps/nwchem/nwchem-7.0.2-release/bin/LINUX64/nwchem'

To Reproduce
Python 3.9.10 | packaged by conda-forge | (main, Jan 30 2022, 18:04:04)
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

import autode as ade
ade.methods.get_lmethod()
MOPAC(available = True)
ade.methods.get_hmethod()
ORCA(available = True)
ade.Config.hcode='NWChem'
ade.methods.get_hmethod()
NWChem(available = True)
h2 = ade.Molecule(smiles='[H][H]')
h2.optimise(method=ade.methods.get_lmethod())
h2.optimise(method=ade.methods.get_hmethod())
Traceback (most recent call last):
File "", line 1, in
File "/home/cnieto/autodE/autode/utils.py", line 227, in wrapped_function
return func(*args, **kwargs)
File "/home/cnieto/autodE/autode/species/species.py", line 826, in optimise
self.atoms = calc.get_final_atoms()
File "/home/cnieto/autodE/autode/calculation.py", line 257, in get_final_atoms
raise ex.AtomsNotFound
autode.exceptions.AtomsNotFound

Expected behavior
reproduce Quick EST Test results (H2)

Environment

• OS: Ubuntu
• Version: 20.04.3 LTS

Additional context

Describe the bug
Attempting to locate a beta-hydride elimination transition state gives "autode.log.log: CRITICAL Reaction had no products – cannot form a reaction"

To Reproduce
Run "python input.py" where input.py is as below:

import autode as ade
r = ade.Reactant(name='reactant',charge=1,mult=4,smiles='CC(C)N1[P+](c2ccccc2)(c3ccccc3)[Cr]4(CCCC4)[P+]1(c5ccccc5)c6ccccc6')
p = ade.Reactant(name='product',charge=1,mult=4,smiles='CC(C)N1[P+](c2ccccc2)(c3ccccc3)[CrH]45(CCC4C5)[P+]1(c6ccccc6)c7ccccc7')
#r.optimise(method=ade.methods.get_lmethod())
#r.energy
#r.atoms
#p.optimise(method=ade.methods.get_lmethod())
#p.energy
#p.atoms
reaction = ade.Reaction(r, p, name='BHE')
print("reaction set up, now calculating reaciton profile")
reaction.calculate_reaction_profile()  # creates BHE/ and saves profile

And I get as output:

autode.log.log: CRITICAL Reaction had no products – cannot form a reaction
Traceback (most recent call last):
  File "autodE_input.py", line 17, in <module>
    reaction = ade.Reaction(r, p, name='BHE')
  File "/home/brc2-skynet/anaconda3/lib/python3.8/site-packages/autode/reactions/reaction.py", line 538, in __init__
    self.type = reaction_types.classify(self.reacs, self.prods)
  File "/home/brc2-skynet/anaconda3/lib/python3.8/site-packages/autode/reactions/reaction_types.py", line 50, in classify
    raise ReactionFormation**Falied**
autode.exceptions.ReactionFormation**Falied**

(minor spelling error bolded above)

If I uncomment the lines for geometry optimization, the corresponding xtb optimizations are performed with no issue, giving reasonable geometries. If I define r and p with the SMILES from the github README, the code works fine:

r = ade.Reactant(name='reactant', smiles='CC[C]([H])[H]')
p = ade.Product(name='product', smiles='C[C]([H])C')

Expected behavior
I expect a transition state search to start.

Environment

• OS: Ubuntu Linux
• Version: 20.04

Additional context

The expected transition state could be quite high in energy, but has been found before for a simplified ligand. See 10.1021/acs.organomet.8b00578 , species "4A-TS3-4" . Thank you for sharing this very clever and intuitive code with the community :)

Describe the bug
Hi there. Is there any way to map atoms implemented currently? When using atom indices as is legal in SMILES and is accepted in Rdkit, I receive the following error:

autode.exceptions.InvalidSmilesString: C[Br:1] had invalid characters:[':']

This would be useful for enforcing/enabling the computation of isomorphic reactions like [Br-:1].C[Br:2]>C[Br:1].[Br-:2]

To Reproduce
molecule = ade.Molecule(smiles='C[Br:1]')

Expected behavior
Generation of molecules with fixed atom indices.

Environment

• Operating System: CentOS-7
• Python version: 3.9.10
• autodE version: 1.2.1

Describe the bug

I want to use the find_tss，but it gives an error,what shou I make my input
To Reproduce

from autode import Reactant, Product, Config, Reaction

from autode.transition_states.locate_tss import find_tss
p = Product('2.xyz')
r = Reactant('1.xyz')
reaction = Reaction(r, p, name='cope')
ts = find_tss(reaction)
Traceback (most recent call last):
File "", line 1, in
File "/data/home/ybw/anaconda3/envs/gsm/lib/python3.9/site-packages/autode/transition_states/locate_tss.py", line 36, in find_tss
if species_are_isomorphic(reactant, product):
File "/data/home/ybw/anaconda3/envs/gsm/lib/python3.9/site-packages/autode/mol_graphs.py", line 219, in species_are_isomorphic
logger.info(f'Checking if {species1.name} and {species2.name} are '
AttributeError: 'NoneType' object has no attribute 'name'

Expected behavior

I want to use the fing_tss
Environment

• OS:Unix
• Version:1.0.2

Additional context

Hi, I would like to run the calculate_energy_profile() method, but it is too expensive (with ORCA) for the reactions I intend to analyze. Since I don't have access to G16 in my machine, I would like to obtain these profiles with XTB as hmethod. However, I've found only the reference for doing it with G16 (with xtb-gaussian):

>>> ade.Config.G16.keywords.sp = ["External='xtb-gaussian'", "IOp(3/5=30)"]
>>> ade.Config.G16.keywords.low_opt = ["External='xtb-gaussian'", "Opt(Loose, NoMicro)", "IOp(3/5=30)"]
>>> ade.Config.G16.keywords.opt = ["External='xtb-gaussian'", "Opt(NoMicro)", "IOp(3/5=30)"]
>>> ade.Config.G16.keywords.opt_ts = ["External='xtb-gaussian'", "Opt(TS, CalcFC, NoEigenTest, MaxCycles=100, MaxStep=10, NoTrustUpdate, NoMicro)", "IOp(3/5=30)"]
>>> ade.Config.G16.keywords.hess = ["External='xtb-gaussian'", "Freq", "Geom(Redundant)", "IOp(3/5=30)"]
>>> ade.Config.G16.keywords.grad = ["External='xtb-gaussian'", 'Force(NoStep)', "IOp(3/5=30)"]

Describe the bug

Some metallocycles e.g. a few from here can raise a ValueError

To Reproduce

import autode as ade
mol = ade.Molecule(smiles='C1[Fe]CCC1')

gives

...autode_env/lib/python3.9/site-packages/autode/smiles/angles.py", line 223, in value
     raise ValueError('Cannot calculate a dihedral - one zero vector')
ValueError: Cannot calculate a dihedral - one zero vector

Expected behavior

These complexes should be build-able!

Environment

• OS: CentOS
• Version: 1.0.3

Additional context

None

Hi! I faced a problem running autodE on my local machine. I run mac and have xtb and Q-Chem installed. I can run calculations with either and they work. Just optimization of a molecule (e.g. water) works well with autodE. However, unfortunately when I try to do "Non-covalent Interaction Complexes" or "Transition States" from your examples, it crashes giving the following error:

multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/opt/anaconda3/envs/dft/lib/python3.9/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/conformers/conformers.py", line 15, in _calc_conformer
    getattr(conformer, calc_type)(method=method,
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/conformers/conformer.py", line 90, in optimise
    super().optimise(method, keywords=keywords, calc=calc, n_cores=n_cores)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/utils.py", line 284, in wrapped_function
    return func(*args, **kwargs)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/species/species.py", line 1043, in optimise
    self.energy = calc.get_energy()
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/calculation.py", line 38, in wrapped_function
    return func(*args, **kwargs)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/calculation.py", line 281, in get_energy
    return PotentialEnergy(self.method.get_energy(self),
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/wrappers/XTB.py", line 199, in get_energy
    raise CouldNotGetProperty(name='energy')
autode.exceptions.CouldNotGetProperty: Could not get energy
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/Users/eremin/Dropbox/Mac/Documents/autodE/autode_test_local7/water_trimer.py", line 23, in <module>
    trimer.find_lowest_energy_conformer(lmethod=xtb)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/utils.py", line 162, in wrapped_function
    result = func(*args, **kwargs)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/species/species.py", line 1213, in find_lowest_energy_conformer
    self.conformers.optimise(method=lmethod)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/conformers/conformers.py", line 256, in optimise
    return self._parallel_calc('optimise', method, keywords)
  File "/opt/anaconda3/envs/dft/lib/python3.9/site-packages/autode/conformers/conformers.py", line 239, in _parallel_calc
    self[idx] = res.get(timeout=None)
  File "/opt/anaconda3/envs/dft/lib/python3.9/multiprocessing/pool.py", line 771, in get
    raise self._value
autode.exceptions.CouldNotGetProperty: Could not get energy

I use the following python script:

import autode as ade
from autode.input_output import xyz_file_to_atoms

xtb = ade.methods.XTB()
qchem = ade.methods.QChem()

ade.Config.n_cores = 4
ade.Config.max_core = 1000

ade.Config.num_complex_sphere_points = 5                                 # N_s

# and the number of rotations to perform per point on the sphere
ade.Config.num_complex_random_rotations = 3                              # N_r

# Total number of conformers ~(N_s × N_r)^(N-1) for N molecules => ~225

# Make a water molecule and optimise at the XTB level
water = ade.Molecule(name='water', smiles='O')
water.optimise(method=xtb)

# Make the NCI complex and find the lowest energy structure
trimer = ade.NCIComplex(water, water, water, name='water_trimer')
trimer.find_lowest_energy_conformer(lmethod=xtb)
trimer.print_xyz_file()

If I run the same script on USC HPC it goes smoothly. This is likely the issue from my system, however, I do not seem to troubleshoot. If someone can advice, it would awesome.
My system:
macOS Monterey 12.3.1
2.9 GHz 6-Core Intel Core i9
32 GB 2400 MHz DDR4
conda 4.12.0
xtb 6.4.1
Q-Chem 5.4.2

Thank you in advance!

Describe the bug
Hi there. In the output folder, file names are given as sum formulas. This creates issues as soon as we look into isomerizations. In the csv file there is no way of telling which of the identically labelled line is what. But much more problematic: both in "reactants_and_products" and in "output" only one set of calculations/one structure will get saved. A quick hotfix would be to check initally if several molecules of the same sum formula exist and add a "b" or a "_1" or something to all filenames.

To Reproduce
This will happen also in the example that is given in the Github manual.

import autode as ade
ade.Config.n_cores = 8

r = ade.Reactant(name='reactant', smiles='CC[C]([H])[H]')
p = ade.Product(name='product', smiles='C[C]([H])C')

reaction = ade.Reaction(r, p, name='1-2_shift')
reaction.calculate_reaction_profile()  # creates 1-2_shift/ and saves profile

Expected behavior
Saving all calculations and structures, ideally in a distinguishable fashion.

Environment

Operating System: CentOS-7
Python version: 3.9.10
autodE version: 1.2.1

Consider dropping support for Python 3.7, which is EOL in 9 months anyway

I have been using the autodE wit G16 and XTB. When I set the following commands, the XTB is executed with 40 cores (2*n_cores) and the Gaussian 16 is executed with 200000 MB (10*max_core) of memory.

Config.n_cores = 20
Config.max_core = 20000

Dear friends ,I have try many reactions, all fails. why the simple CH4+S--CSH2+H2 can not rightly local?

s-h2.zip

Describe the bug

XTB calculations are not run using the spin multiplicity defined in the molecule. Missing "--uhf N" where N is the number of unpaired electrons

To Reproduce

>>> import autode as ade
>>> mol = ade.Molecule(smiles='C[C]C', mult=3)
>>> mol.mult
3
>>> mol.single_point(method=ade.methods.XTB())

the XTB output file contains:

*****
         program call                      : /../bin/xtb C3H6_sp_xtb.xyz --chrg 0
                                ****
          spin                             : 0.0
*****

Expected behavior

XTB should be called with the additional --uhf mol.mult-1 flag

Environment

• OS: Linux
• Version: 1.0.5

Hello again,

just a tiny question as I am putting together the SI for my first paper utilising autodE and I have encountered a little question. I have some difficulties finding out what exactly happens to imaginary frequencies that are cut off due to being lower than -40 cm**-1 in the thermochemical calculations. Are they treated as positive or as zero?

Thank you!

Hello and thank you for the package.
I'm not super knowledgeable about python or conda, but for a quickstart installation guide, maybe the following is useful for others:

I found nwchem as the easiest install for the first dependency requirement, ORCA seemed to require making an account somewhere.

conda create -n myautode python=3.9
conda install autode --channel conda-forge
conda install ipykernel
conda install xtb
# resolution issue if using python 3.10 as of 5/30/22
conda install -c conda-forge nwchem

looking forward to using autode! feel free to close the issue

Describe the bug
I´d like to restart a TS search of a reaction invoked as rxn.calculate_reaction_profile(). It stopped for some memory issues (not autode related) and I´d like to avoid the hassle of recalculating conformers and reactants and products. Is there any specific commands to do this?

Thank you kindly.

E77

To Reproduce

Expected behavior

Environment

• Operating System:
• Python version:
• autodE version:

Additional context

Describe the bug
Adding solvent variables seems to function differently in transition states than in other parts of the code. I understand that this usually will take its variables from the provided reactant and product but it is possible to use the script to run a standard TS optimisation straight from an xyz-file without reactant and product complexes. This is helpful for enabling workarounds for reactions that autodE cannot do at the moment (i.e. I am doing the rotation barrier calculations that I mentioned previously via manually guessed structures or via PES-scan guesses). Neither TSguess() nor TransitionState() eats "solvent_name". TSguess() does take, however, "charge" and "mult" so it appears inconsistent with how Molecule() works that "solvent_name" cannot be used. TransitionState() does take TransitionState().solvent, though (which is then the only part that was defined on this function and not on TSguess).

To Reproduce
Assume we have a molecule defined.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult,
                                   solvent_name=molecule.solvent
                                   )
molecule = TransitionState(tsguess)
molecule.optimise(method=hmethod)

This throws an unexpected keyword.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
molecule = TransitionState(tsguess, solvent_name = molecule.solvent)
molecule.optimise(method=hmethod)

This throws an unexpected keyword.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
molecule = TransitionState(tsguess)
molecule.solvent = get_solvent(molecule.solvent, kind='implicit')
molecule.optimise(method=hmethod)

This works.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
tsguess.solvent = get_solvent(molecule.solvent, kind='implicit')
molecule = TransitionState(tsguess)
molecule.optimise(method=hmethod)

This one is dangerous, it did not throw an exception but the CPCM(acetonitrile) is nowhere to be found in the generated ORCA inputs. When my reaction energy was roughly 700 kJ/mol wrong, I got a little suspicious here.

Expected behavior
Solvents should be possible to define where charge and mult is allowed as well.