diff --git a/carmm/analyse/counterpoise_onepot.py b/carmm/analyse/counterpoise_onepot.py index 0e5a94c0..5f9fafd9 100644 --- a/carmm/analyse/counterpoise_onepot.py +++ b/carmm/analyse/counterpoise_onepot.py @@ -57,6 +57,11 @@ def counterpoise_calc(complex_struc, a_id, b_id, fhi_calc=None, a_name=None, b_n # Empty sites does not work with forces. Remove compute_forces. if 'compute_forces' in fhi_calc.parameters: fhi_calc.parameters.pop('compute_forces') + if 'sc_accuracy_forces' in fhi_calc.parameters: + if verbose: + print('Stop calculation as there is a convergence criterion regarding force.', '\n', + 'Empty sites does not work with forces. Remove and check how it affects your results.') + raise KeyError('Found sc_accuracy_forces in parameters! FHI-aims can not calculate forces on empty sites.') # Create an empty list to store energies for postprocessing. energies = [] for index in range(4): @@ -72,10 +77,9 @@ def counterpoise_calc(complex_struc, a_id, b_id, fhi_calc=None, a_name=None, b_n energies.append(energy_i) else: fhi_calc.template.outputname = species_list[index] + '.out' - properties = fhi_calc.implemented_properties - parameters = fhi_calc.parameters - parameters['ghosts'] = ghosts_lists_cp[index] - fhi_calc.template.update_parameters(properties, parameters) + fhi_calc.parameters['ghosts'] = ghosts_lists_cp[index] + # Scaled positions does not work with empty sites. + fhi_calc.parameters['scaled'] = False structures_cp[index].calc = fhi_calc if dry_run: structures_cp[index].calc.template.write_input(fhi_calc.profile, fhi_calc.directory, @@ -192,14 +196,17 @@ def calculate_energy_ghost_compatible(calc, atoms=None, properties=['energy'], """ from ase.calculators.calculator import Calculator - import subprocess + import subprocess, os Calculator.calculate(calc, atoms, properties, system_changes) - calc.write_input(calc.atoms, properties, system_changes, ghosts=ghosts) + # Write inputfiles. Scaled positions does not work with empty sites. + calc.write_input(calc.atoms, properties, system_changes, ghosts=ghosts, scaled=False) command = calc.command if dry_run: # Only for CI tests command = '' # Used to be 'ls' - converged = calc.read_convergence() + converged = False + if os.path.exists(calc.directory+'/'+calc.outfilename): + converged = calc.read_convergence() if (not converged) or dry_run: subprocess.check_call(command, shell=True, cwd=calc.directory) diff --git a/examples/analyse_counterpoise.py b/examples/analyse_counterpoise.py index 93faf352..f089d7a6 100644 --- a/examples/analyse_counterpoise.py +++ b/examples/analyse_counterpoise.py @@ -14,6 +14,7 @@ def test_analyse_counterpoise(): from carmm.analyse.counterpoise_onepot import counterpoise_calc from carmm.run.aims_calculator import get_aims_calculator from carmm.run.aims_path import set_aims_command + from unittest import TestCase # This is an example script for using counterpoise_calc for counterpoise (CP) correction. Please note the species # files in data/CO_BSSE are fake ones and default species settings are also deleted from aims.out. @@ -39,6 +40,12 @@ def test_analyse_counterpoise(): verbose=True, dry_run=True) cp_symbol = counterpoise_calc(CO, a_id=['C'], b_id=['O'], fhi_calc=toy_calc, dry_run=True) + # Test to make sure sc_accuracy_forces is not set, as FHI-aims can't calculate forces with empty sites + toy_calc.parameters['sc_accuracy_forces'] = 0.0001 + forces_check = TestCase() + with forces_check.assertRaises(KeyError): + counterpoise_calc(CO, a_id=['C'], b_id=['O'], fhi_calc=toy_calc, dry_run=True, verbose=True) + # CP correction = A_only + B_only - A_plus_ghost - B_plus_ghost # This value should be added to the energy change of interest, such as adsorption energy.