Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changed the uhf value when xtb calculates with lanthanides #60

Merged
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

### Breaking Changes
- Removal of the `dist_threshold` flag and in the `-toml` file.
- The number of unpaired electrons (`Molecule.uhf`) is now set to 0 if `xtb` is used as `QMMethod` and a lanthanide is within the molecule to match the `f-in-core` approximation.

## [0.4.0] - 2024-09-19
### Changed
Expand Down
45 changes: 43 additions & 2 deletions src/mindlessgen/qm/xtb.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,11 @@
from pathlib import Path
import shutil
from tempfile import TemporaryDirectory
import numpy as np
from ..molecules import Molecule
from .base import QMMethod
from ..prog import XTBConfig
from ..molecules.miscellaneous import get_lanthanides


class XTB(QMMethod):
Expand All @@ -34,7 +36,24 @@ def optimize(
"""
Optimize a molecule using xtb.
"""

if np.any(np.isin(molecule.ati, get_lanthanides())):
nel = 0
f_electrons = 0
for atom in molecule.ati:
nel += atom + 1
if atom in get_lanthanides():
f_electrons += atom - 56
# Check if the number of the remaning electrons is even
test_rhf = nel - f_electrons - molecule.charge
if not test_rhf % 2 == 0:
raise SystemError(
"Lanthanides detected, remaining number of electrons is not even."
)
# Store the original UHF value and set uhf to 0
# Justification: xTB does not treat f electrons explicitly.
# The remaining openshell system has to be removed.
uhf_original = molecule.uhf
molecule.uhf = 0
# Create a unique temporary directory using TemporaryDirectory context manager
with TemporaryDirectory(prefix="xtb_") as temp_dir:
temp_path = Path(temp_dir).resolve()
Expand Down Expand Up @@ -71,13 +90,33 @@ def optimize(
# read the optimized molecule
optimized_molecule = molecule.copy()
optimized_molecule.read_xyz_from_file(temp_path / "xtbopt.xyz")

if np.any(np.isin(molecule.ati, get_lanthanides())):
# Reset the UHF value to the original value before returning the optimized molecule.
optimized_molecule.uhf = uhf_original
return optimized_molecule

def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str:
"""
Perform a single-point calculation using xtb.
"""
if np.any(np.isin(molecule.ati, get_lanthanides())):
nel = 0
f_electrons = 0
for atom in molecule.ati:
nel += atom + 1
if atom in get_lanthanides():
f_electrons += atom - 56
# Check if the number of the remaning electrons is even
test_rhf = nel - f_electrons - molecule.charge
if not test_rhf % 2 == 0:
raise SystemError(
"Lanthanides detected, remaining number of electrons is not even."
)
jonathan-schoeps marked this conversation as resolved.
Show resolved Hide resolved
# Store the original UHF value and set uhf to 0
# Justification: xTB does not treat f electrons explicitly.
# The remaining openshell system has to be removed.
uhf_original = molecule.uhf
molecule.uhf = 0

# Create a unique temporary directory using TemporaryDirectory context manager
with TemporaryDirectory(prefix="xtb_") as temp_dir:
Expand Down Expand Up @@ -109,6 +148,8 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str:
f"xtb failed with return code {return_code}:\n{xtb_log_err}"
)

if np.any(np.isin(molecule.ati, get_lanthanides())):
molecule.uhf = uhf_original
return xtb_log_out

def check_gap(
Expand Down