From 660206521bb2ff07cde27d86b5360119ecd59a5c Mon Sep 17 00:00:00 2001 From: tom Date: Tue, 31 Jan 2023 21:10:57 +0000 Subject: [PATCH] Fix tests --- autode/neb/original.py | 5 ++-- autode/opt/optimisers/dimer.py | 5 ++-- autode/transition_states/locate_tss.py | 3 +-- doc/common/nci_FF_example.py | 14 +++++------ tests/test_path.py | 2 ++ tests/test_pes/test_calculate.py | 2 +- tests/test_ts/test_mode_checking.py | 7 +++--- tests/test_ts/test_ts_adapt_neb.py | 3 ++- tests/test_wrappers/test_gaussian.py | 26 ++++++++++++--------- tests/test_wrappers/test_mopac.py | 32 +++++++++++--------------- tests/test_wrappers/test_orca.py | 26 ++++++++++----------- tests/test_wrappers/test_qchem.py | 24 ++++++++----------- 12 files changed, 73 insertions(+), 76 deletions(-) diff --git a/autode/neb/original.py b/autode/neb/original.py index 1de82d1c9..1787f54c7 100644 --- a/autode/neb/original.py +++ b/autode/neb/original.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib.pyplot as plt -from typing import Optional, Sequence, List, Any, Union, TYPE_CHECKING +from typing import Optional, Sequence, List, Any, TYPE_CHECKING from copy import deepcopy from autode.log import logger @@ -168,7 +168,8 @@ def __init__( mult=species.mult, atoms=species.atoms.copy(), ) - self.solvent = species.solvent + self.solvent = deepcopy(species.solvent) + self.energy = deepcopy(species.energy) self.iteration = 0 #: Current optimisation iteration of this image self.k = k diff --git a/autode/opt/optimisers/dimer.py b/autode/opt/optimisers/dimer.py index 92e592b17..cefd8f15c 100644 --- a/autode/opt/optimisers/dimer.py +++ b/autode/opt/optimisers/dimer.py @@ -273,11 +273,10 @@ def _update_gradient_at(self, point: DimerPoint) -> None: ) calc.run() - self._coords.e = self._species.energy = calc.get_energy() - self._species.gradient = calc.get_gradients() + self._coords.e = self._species.energy self._coords.set_g_at( - point, calc.get_gradients().flatten(), mass_weighted=False + point, self._species.gradient.flatten(), mass_weighted=False ) calc.clean_up(force=True, everything=True) diff --git a/autode/transition_states/locate_tss.py b/autode/transition_states/locate_tss.py index 00d1c5fb2..e1163f834 100644 --- a/autode/transition_states/locate_tss.py +++ b/autode/transition_states/locate_tss.py @@ -356,9 +356,8 @@ def _get_ts_neb_from_adaptive_path( ) if neb.images.contains_peak: - peak_idx = neb.images.peak_idx ts_guess = TSguess( - atoms=neb.images[peak_idx].species.atoms, + atoms=neb.peak_species.atoms, reactant=reactant, product=product, bond_rearr=bond_rearr, diff --git a/doc/common/nci_FF_example.py b/doc/common/nci_FF_example.py index e966b30ba..0d865d652 100644 --- a/doc/common/nci_FF_example.py +++ b/doc/common/nci_FF_example.py @@ -1,8 +1,7 @@ import autode as ade -from autode.methods import XTB -from autode.calculations import Calculation -from scipy.optimize import minimize import numpy as np + +from scipy.optimize import minimize from scipy.spatial import distance_matrix ade.Config.max_num_complex_conformers = 10 @@ -85,17 +84,16 @@ def rotation_translation( def set_charges_vdw(species): """Calculate the partial atomic charges to atoms with XTB""" - calc = Calculation( + calc = ade.Calculation( name="tmp", molecule=species, - method=XTB(), - keywords=ade.SinglePointKeywords([]), + method=ade.methods.XTB(), + keywords=ade.SinglePointKeywords(), ) calc.run() - charges = calc.get_atomic_charges() for i, atom in enumerate(species.atoms): - atom.charge = charges[i] + atom.charge = atom.partial_charge atom.vdw = float(atom.vdw_radius) return None diff --git a/tests/test_path.py b/tests/test_path.py index 32a62abd2..237beffa6 100644 --- a/tests/test_path.py +++ b/tests/test_path.py @@ -9,6 +9,7 @@ from autode.bonds import FormingBond, BreakingBond from autode.species import Species, Molecule from autode.units import Unit, KcalMol +from . import testutils test_species = Species(name="tmp", charge=0, mult=1, atoms=[Atom("He")]) test_mol = Molecule(smiles="O") @@ -179,6 +180,7 @@ def test_products_made(): assert not path.products_made(product=diff_mol) +@testutils.requires_with_working_xtb_install def test_adaptive_path(): species_no_atoms = Species(name="tmp", charge=0, mult=1, atoms=[]) diff --git a/tests/test_pes/test_calculate.py b/tests/test_pes/test_calculate.py index f2abf9bf0..3f2f3ee4d 100644 --- a/tests/test_pes/test_calculate.py +++ b/tests/test_pes/test_calculate.py @@ -24,7 +24,7 @@ def test_calculate_no_species(): @testutils.requires_with_working_xtb_install @work_in_tmp_dir(filenames_to_copy=[], kept_file_exts=[]) def test_calculate_1d(): - pes = RelaxedPESnD(species=h2(), rs={(0, 1): (1.6, 10)}) + pes = RelaxedPESnD(species=h2(), rs={(0, 1): (1.5, 10)}) pes.calculate(method=XTB()) diff --git a/tests/test_ts/test_mode_checking.py b/tests/test_ts/test_mode_checking.py index 702c9f984..199386287 100644 --- a/tests/test_ts/test_mode_checking.py +++ b/tests/test_ts/test_mode_checking.py @@ -81,7 +81,8 @@ def has_correct_mode(name, fbonds, bbonds): keywords=orca.keywords.opt_ts, n_cores=1, ) - calc.molecule = reac = Reactant( + # need to bypass the pre-calculation checks on the molecule. e.g. valid spin state + calc.molecule = reactant = Reactant( name="r", atoms=xyz_file_to_atoms(f"{name}.xyz") ) @@ -89,11 +90,11 @@ def has_correct_mode(name, fbonds, bbonds): # Don't require all bonds to be breaking/making in a 'could be ts' function ts = TSbase( - atoms=reac.atoms, + atoms=reactant.atoms, bond_rearr=BondRearrangement( breaking_bonds=bbonds, forming_bonds=fbonds ), ) - ts.hessian = calc.get_hessian() + ts.hessian = reactant.hessian return ts.imag_mode_has_correct_displacement(req_all=False) diff --git a/tests/test_ts/test_ts_adapt_neb.py b/tests/test_ts/test_ts_adapt_neb.py index bff7319ab..448d92edb 100644 --- a/tests/test_ts/test_ts_adapt_neb.py +++ b/tests/test_ts/test_ts_adapt_neb.py @@ -43,8 +43,9 @@ def test_ts_from_neb_optimised_after_adapt(): rxn = _sn2_reaction() neb = NEB.from_file("dOr2Us_ll_ad_0-5_0-1_path.xyz") + assert len(neb.images) > 0 - assert neb.images[0].species.solvent is not None + assert neb.images[0].solvent is not None def get_ts_guess(): return _get_ts_neb_from_adaptive_path( diff --git a/tests/test_wrappers/test_gaussian.py b/tests/test_wrappers/test_gaussian.py index cd8310e54..3b96756de 100644 --- a/tests/test_wrappers/test_gaussian.py +++ b/tests/test_wrappers/test_gaussian.py @@ -23,7 +23,6 @@ from .. import testutils here = os.path.dirname(os.path.abspath(__file__)) -test_mol = Molecule(name="methane", smiles="C") method = G09() opt_keywords = OptKeywords(["PBE1PBE/Def2SVP", "Opt"]) @@ -39,6 +38,10 @@ sp_keywords = SinglePointKeywords(["PBE1PBE/Def2SVP"]) +def methane(): + return Molecule(name="methane", smiles="C") + + def test_printing_ecp(): tmp_file = open("tmp.com", "w") @@ -90,7 +93,7 @@ def test_input_print_max_opt(): keywds = opt_keywords.copy() keywds.max_opt_cycles = 10 - str_keywords = _get_keywords(CalculationInput(keywds), molecule=test_mol) + str_keywords = _get_keywords(CalculationInput(keywds), molecule=methane()) # Should be only a single instance of the maxcycles declaration assert sum("maxcycles=10" in kw.lower() for kw in str_keywords) == 1 @@ -215,7 +218,7 @@ def test_bad_gauss_output(): calc = Calculation( name="no_output", - molecule=test_mol, + molecule=methane(), method=method, keywords=opt_keywords, ) @@ -248,7 +251,7 @@ def test_fix_angle_error(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "g09.zip")) def test_constraints(): - a = test_mol.copy() + a = methane() a.constraints.distance = {(0, 1): 1.2} calc = Calculation( name="const_dist_opt", molecule=a, method=method, keywords=opt_keywords @@ -260,14 +263,14 @@ def test_constraints(): 1.199 < np.linalg.norm(opt_atoms[0].coord - opt_atoms[1].coord) < 1.201 ) - b = test_mol.copy() + b = methane() b.constraints.cartesian = [0] calc = Calculation( name="const_cart_opt", molecule=b, method=method, keywords=opt_keywords ) calc.run() opt_atoms = b.atoms - assert np.linalg.norm(test_mol.atoms[0].coord - opt_atoms[0].coord) < 1e-3 + assert np.linalg.norm(methane().atoms[0].coord - opt_atoms[0].coord) < 1e-3 @testutils.work_in_zipped_dir(os.path.join(here, "data", "g09.zip")) @@ -299,9 +302,10 @@ def test_point_charge_calc(): # Methane single point using a point charge with a unit positive charge # located at (10, 10, 10) + mol = methane() calc = Calculation( name="methane_point_charge", - molecule=test_mol, + molecule=mol, method=method, keywords=sp_keywords, point_charges=[PointCharge(charge=1.0, x=10.0, y=10.0, z=10.0)], @@ -323,14 +327,14 @@ def test_point_charge_calc(): assert float(z) == 10.0 assert float(charge) == 1.0 - assert -40.428 < calc.get_energy() < -40.427 + assert -40.428 < mol.energy < -40.427 # Gaussian needs x-matrix and nosymm in the input line to run optimisations # with point charges.. for opt_keyword in ["Opt", "Opt=Tight", "Opt=(Tight)"]: calc = Calculation( name="methane_point_charge_o", - molecule=test_mol, + molecule=methane(), method=method, keywords=OptKeywords(["PBE1PBE/Def2SVP", opt_keyword]), point_charges=[PointCharge(charge=1.0, x=3.0, y=3.0, z=3.0)], @@ -417,5 +421,5 @@ def test_xtb_optts(): calc.run() # Even though a Hessian is not requested it should be added - assert calc.get_hessian() is not None - assert np.isclose(calc.get_energy().to("Ha"), -13.1297380, atol=1e-5) + assert orca_ts.hessian is not None + assert np.isclose(orca_ts.energy.to("Ha"), -13.1297380, atol=1e-5) diff --git a/tests/test_wrappers/test_mopac.py b/tests/test_wrappers/test_mopac.py index e137c99b8..e4cc07d68 100644 --- a/tests/test_wrappers/test_mopac.py +++ b/tests/test_wrappers/test_mopac.py @@ -92,7 +92,6 @@ def test_other_spin_states(): keywords=Config.MOPAC.keywords.sp, ) calc.run() - singlet_energy = calc.get_energy() o_triplet = Molecule(atoms=[Atom("O")], mult=3) o_triplet.name = "molecule" @@ -104,9 +103,8 @@ def test_other_spin_states(): keywords=Config.MOPAC.keywords.sp, ) calc.run() - triplet_energy = calc.get_energy() - assert triplet_energy < singlet_energy + assert o_triplet.energy < o_singlet.energy h_doublet = Molecule(atoms=[Atom("H")], mult=2) h_doublet.name = "molecule" @@ -120,7 +118,7 @@ def test_other_spin_states(): calc.run() # Open shell doublet should work - assert calc.get_energy() is not None + assert h_doublet.energy is not None h_quin = Molecule(atoms=[Atom("H")], mult=5) h_quin.name = "molecule" @@ -148,11 +146,9 @@ def test_bad_geometry(): keywords=Config.MOPAC.keywords.opt, ) - calc.output.filename = "h2_overlap_opt_mopac.out" - assert not calc.terminated_normally - with pytest.raises(Exception): - _ = calc.get_energy() + # cannot even get the energy from the output file + calc.set_output_filename("h2_overlap_opt_mopac.out") assert not method.optimiser_from(calc).converged @@ -169,21 +165,23 @@ def test_constrained_opt(): keywords=Config.MOPAC.keywords.opt, ) calc.run() - opt_energy = calc.get_energy() + opt_energy = methane.energy # Constrained optimisation with a C–H distance of 1.2 Å # (carbon is the first atom in the file) + constrained_methane = methane.copy() methane.constraints.distance = {(0, 1): 1.2} + const = Calculation( name="methane_const", - molecule=methane, + molecule=constrained_methane, method=method, keywords=Config.MOPAC.keywords.opt, ) const.run() - assert opt_energy < const.get_energy() - assert calc.get_hessian() is None + assert methane.energy < constrained_methane.energy + assert methane.hessian is None @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) @@ -198,10 +196,9 @@ def test_grad(): keywords=Config.MOPAC.keywords.grad, ) grad_calc.run() - energy = grad_calc.get_energy() - assert energy is not None + assert h2.energy is not None - gradients = grad_calc.get_gradients() + gradients = h2.gradient assert gradients.shape == (2, 3) delta_r = 1e-5 @@ -210,7 +207,7 @@ def test_grad(): ) h2_disp.single_point(method) - delta_energy = h2_disp.energy - energy # Ha] + delta_energy = h2_disp.energy - h2.energy grad = delta_energy / delta_r # Ha A^-1 # Difference between the absolute and finite difference approximation @@ -227,10 +224,9 @@ def test_broken_grad(): method=method, keywords=Config.MOPAC.keywords.grad, ) - grad_calc_broken.output.filename = "h2_grad_broken.out" with pytest.raises(CouldNotGetProperty): - _ = grad_calc_broken.get_gradients() + grad_calc_broken.set_output_filename("h2_grad_broken.out") @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) diff --git a/tests/test_wrappers/test_orca.py b/tests/test_wrappers/test_orca.py index 1f2b78a89..2b5734ccb 100644 --- a/tests/test_wrappers/test_orca.py +++ b/tests/test_wrappers/test_orca.py @@ -238,7 +238,7 @@ def test_gradients(): calc.run() - diff = calc.get_gradients()[1, 0] - grad # Ha A^-1 + diff = h2.gradient[1, 0] - grad # Ha A^-1 # Difference between the absolute and finite difference approximation assert np.abs(diff) < 1e-3 @@ -247,15 +247,16 @@ def test_gradients(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "orca.zip")) def test_mp2_numerical_gradients(): + mol = Molecule("tmp_orca.xyz", charge=-1) calc = Calculation( name="tmp", - molecule=Molecule(atoms=xyz_file_to_atoms("tmp_orca.xyz"), charge=-1), + molecule=mol, method=method, keywords=method.keywords.grad, ) calc.set_output_filename(filename="tmp_orca.out") - gradients = calc.get_gradients() + gradients = mol.gradient assert len(gradients) == 6 expected = ( np.array([-0.00971201, -0.00773534, -0.02473580]) / Constants.a0_to_ang @@ -266,7 +267,7 @@ def test_mp2_numerical_gradients(): calc.set_output_filename(filename="numerical_orca.out") assert calc.output.filename == "numerical_orca.out" - gradients = calc.get_gradients() + gradients = mol.gradient assert len(gradients) == 6 expected = ( np.array([0.012397372, 0.071726232, -0.070942743]) @@ -323,15 +324,14 @@ def test_keyword_setting(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "orca.zip")) def test_hessian_extraction(): + h2o = Molecule(smiles="O") calc = Calculation( name="tmp", - molecule=Molecule(smiles="O"), + molecule=h2o, method=method, keywords=method.keywords.hess, ) - calc.output.filename = "H2O_hess_orca.out" - with open("H2O_hess_orca.xyz", "w") as xyz_file: print( "3\n", @@ -342,7 +342,9 @@ def test_hessian_extraction(): file=xyz_file, ) - hessian = calc.get_hessian() + calc.set_output_filename("H2O_hess_orca.out") + + hessian = h2o.hessian assert hessian.shape == (9, 9) # should not have any very large values assert np.sum(np.abs(hessian)) < 100 @@ -364,12 +366,10 @@ def test_charges_from_v5_output_file(): method=method, keywords=method.keywords.sp, ) - calc.output.filename = "h2o_orca_v5_charges.out" - + calc.set_output_filename("h2o_orca_v5_charges.out") assert calc.output.exists - - # q_O q_H q_H - assert calc.get_atomic_charges() == [-0.313189, 0.156594, 0.156594] + # q_O q_H q_H + assert water.partial_charges == [-0.313189, 0.156594, 0.156594] def test_unsupported_freq_scaling(): diff --git a/tests/test_wrappers/test_qchem.py b/tests/test_wrappers/test_qchem.py index ee842e63d..176c7b1bf 100644 --- a/tests/test_wrappers/test_qchem.py +++ b/tests/test_wrappers/test_qchem.py @@ -35,7 +35,7 @@ def _blank_calc(name="test"): def _completed_thf_calc(): calc = _blank_calc() - calc.output.filename = "smd_thf.out" + calc.set_output_filename("smd_thf.out") assert calc.output.exists assert len(calc.output.file_lines) > 0 @@ -204,9 +204,8 @@ def test_simple_input_generation(): def test_energy_extraction(): calc = _completed_thf_calc() - energy = calc.get_energy() - assert np.isclose(energy.to("Ha"), -232.45463628, atol=1e-8) + assert np.isclose(calc.molecule.energy.to("Ha"), -232.45463628, atol=1e-8) for calc in (_blank_calc(), _broken_output_calc(), _broken_output_calc2()): with pytest.raises(CalculationException): @@ -294,9 +293,11 @@ def test_h2o_opt(): @work_in_zipped_dir(qchem_data_zip_path) def test_gradient_extraction_h2o(): + + h2o = Molecule(smiles="O") calc = Calculation( name="test", - molecule=Molecule(smiles="O"), + molecule=h2o, method=method, keywords=OptKeywords(), ) @@ -305,15 +306,13 @@ def test_gradient_extraction_h2o(): assert calc.output.exists - grad = calc.get_gradients() - assert grad.shape == (3, 3) + assert h2o.gradient.shape == (3, 3) # The minimum should have a gradient close to zero - assert np.allclose(grad, np.zeros(shape=(3, 3)), atol=1e-4) + assert np.allclose(h2o.gradient, np.zeros(shape=(3, 3)), atol=1e-4) # also for this calculation the optimisation has converged assert calc.optimiser.converged - assert not calc.optimisation_nearly_converged() @work_in_zipped_dir(qchem_data_zip_path) @@ -321,10 +320,9 @@ def test_gradient_extraction_h2(): calc = _blank_calc() calc.molecule = Molecule(atoms=[Atom("H"), Atom("H", x=0.77)]) - calc.output.filename = "H2_qchem.out" + calc.set_output_filename("H2_qchem.out") - grad = calc.get_gradients() - assert grad.shape == (2, 3) + assert calc.molecule.gradient.shape == (2, 3) @work_in_zipped_dir(qchem_data_zip_path) @@ -447,9 +445,7 @@ def test_ts_opt(): # Should skip calculation for already completed and saved calculation calc.run() - assert np.isclose(calc.get_energy().to("Ha"), -599.4788133790, atol=1e-8) - - ts_mol.hessian = calc.get_hessian() + assert np.isclose(ts_mol.energy.to("Ha"), -599.4788133790, atol=1e-8) assert ts_mol.hessian is not None assert sum(freq.is_imaginary for freq in ts_mol.vib_frequencies) == 1