diff --git a/CHANGELOG.md b/CHANGELOG.md index a00ab05..098031b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Support for the novel "g-xTB" method (working title: GP3-xTB) +- A function which contracts the coordinates after the initial generation. +- A function which is able to printout the xyz coordinates to the terminal similar to the `.xyz` layout. ### Breaking Changes - Removal of the `dist_threshold` flag and in the `-toml` file. diff --git a/mindlessgen.toml b/mindlessgen.toml index d6501e3..9cb79b5 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -28,9 +28,11 @@ init_scaling = 3.0 # > Increase in the coordinate scaling factor per trial after check_distance was not met. Options: increase_scaling_factor = 1.3 # > Scaling factor for the employed van der Waals radii. Options: -scale_vdw_radii = 1.3333 +scale_vdw_radii = 1.25 # > Scaling factor for the minimal bondlength based on the sum of the van der Waals radii. Options: -scale_minimal_bondlength = 0.75 +scale_minimal_bondlength = 0.8 +# > Contract the coordinates after the initial generation. Leads to more cluster-like and less extended structures. Options: +contract_coords = false # > Atom types and their minimum and maximum occurrences. Format: ":-" # > Elements that are not specified are only added by random selection. # > A star sign (*) can be used as a wildcard for integer value. diff --git a/src/mindlessgen/cli/cli_parser.py b/src/mindlessgen/cli/cli_parser.py index 2a8e345..8480108 100644 --- a/src/mindlessgen/cli/cli_parser.py +++ b/src/mindlessgen/cli/cli_parser.py @@ -137,6 +137,12 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: required=False, help="List of forbidden elements.", ) + parser.add_argument( + "--contract_coords", + type=bool, + required=False, + help="Contract the coordinates of the molecule after the coordinats generation.", + ) ### Refinement arguments ### parser.add_argument( @@ -195,6 +201,12 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: required=False, help="Path to the xTB binary.", ) + parser.add_argument( + "--xtb-level", + type=int, + required=False, + help="Level of theory to use in xTB.", + ) parser.add_argument( "--postprocess-debug", action="store_true", @@ -202,12 +214,6 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: required=False, help="Print debug information during postprocessing.", ) - parser.add_argument( - "--xtb-level", - type=int, - required=False, - help="Level of theory to use in xTB.", - ) ### ORCA specific arguments ### parser.add_argument( @@ -273,6 +279,7 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: "forbidden_elements": args_dict["forbidden_elements"], "scale_vdw_radii": args_dict["scale_vdw_radii"], "scale_minimal_bondlength": args_dict["scale_minimal_bondlength"], + "contract_coords": args_dict["contract_coords"], } # XTB specific arguments rev_args_dict["xtb"] = {"xtb_path": args_dict["xtb_path"]} diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 0f0eac4..9de4d88 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -41,6 +41,12 @@ def generate_random_molecule( verbosity=verbosity, scale_bondlength=config_generate.scale_minimal_bondlength, ) + if config_generate.contract_coords: + mol.xyz = contract_coordinates( + xyz=mol.xyz, + ati=mol.ati, + scale_minimal_distance=config_generate.scale_minimal_bondlength, + ) mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) mol.set_name_from_formula() @@ -371,6 +377,39 @@ def generate_random_coordinates(at: np.ndarray) -> tuple[np.ndarray, np.ndarray] return xyz, ati +def contract_coordinates( + xyz: np.ndarray, ati: np.ndarray, scale_minimal_distance: float +) -> np.ndarray: + """ + Pull the atoms towards the origin. + """ + # Initialize the old coordinates as an array of zeros + xyz_old: np.ndarray = np.zeros_like(xyz) + cycle = 0 + # Break if the coordinates do not change + while not np.array_equal(xyz_old, xyz): + cycle += 1 + if cycle > 2500: + raise RuntimeError( + "Could not contract the coordinates in a reasonable amount of cycles." + ) + xyz_old = xyz.copy() + # Go through the atoms dimension of the xyz array in a reversed order + # Justification: First atoms are most likely hydrogen atoms, which should be moved last + for i in range(len(xyz) - 1, -1, -1): + atom_xyz = xyz[i] + atom_xyz_norm = np.linalg.norm(atom_xyz) + normalized_atom_xyz = atom_xyz / atom_xyz_norm + # Shift the atom only if it is closer to the origin after the shift + if atom_xyz_norm > 0.1: + # Pull the atom towards the origin + xyz[i] -= normalized_atom_xyz * 0.2 + # When the check_distances function returns False, reset the atom coordinates + if not check_distances(xyz, ati, scale_minimal_distance): + xyz[i] = xyz_old[i] + return xyz + + def check_distances(xyz: np.ndarray, ati: np.ndarray, scale_bondlength: float) -> bool: """ Check if the distances between atoms are larger than a threshold. diff --git a/src/mindlessgen/molecules/molecule.py b/src/mindlessgen/molecules/molecule.py index e14f418..33986f8 100644 --- a/src/mindlessgen/molecules/molecule.py +++ b/src/mindlessgen/molecules/molecule.py @@ -459,12 +459,29 @@ def atlist(self, value: np.ndarray): self._atlist = value - def print_xyz(self): + def get_xyz_str(self) -> str: """ - Print the XYZ coordinates of the molecule. + Obtain a string with the full XYZ file information of the molecule. """ - # TODO: Implement a better way to print the coordinates - print(self._xyz) + xyz_str = f"{self.num_atoms}\n" + try: + commentline = f"Total charge: {self.charge} ; " + except ValueError: + commentline = "" + try: + commentline = commentline + f"Unpaired electrons: {self.uhf} ; " + except ValueError: + pass + commentline = commentline + f"Generated by mindlessgen-v{__version__}\n" + xyz_str += commentline + for i in range(self.num_atoms): + xyz_str += ( + f"{PSE[self.ati[i]+1]:<5} " + + f"{self.xyz[i, 0]:>12.7f} " + + f"{self.xyz[i, 1]:>12.7f} " + + f"{self.xyz[i, 2]:>12.7f}\n" + ) + return xyz_str def write_xyz_to_file(self, filename: str | Path | None = None): """ @@ -496,24 +513,7 @@ def write_xyz_to_file(self, filename: str | Path | None = None): filename = Path("mlm_" + self.name + ".xyz").resolve() with open(filename, "w", encoding="utf8") as f: - f.write(f"{self.num_atoms}\n") - try: - commentline = f"Total charge: {self.charge} ; " - except ValueError: - commentline = "" - try: - commentline = commentline + f"Unpaired electrons: {self.uhf} ; " - except ValueError: - pass - commentline = commentline + f"Generated by mindlessgen-v{__version__}\n" - f.write(commentline) - for i in range(self.num_atoms): - f.write( - f"{PSE[self.ati[i]+1]:<5} " - + f"{self.xyz[i, 0]:>12.7f} " - + f"{self.xyz[i, 1]:>12.7f} " - + f"{self.xyz[i, 2]:>12.7f}\n" - ) + f.write(self.get_xyz_str()) # if the charge is set, write it to a '.CHRG' file if self._charge is not None: with open(filename.with_suffix(".CHRG"), "w", encoding="utf8") as f: diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 96ffe68..98688ab 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -178,8 +178,9 @@ def __init__(self: GenerateConfig) -> None: self._increase_scaling_factor: float = 1.3 self._element_composition: dict[int, tuple[int | None, int | None]] = {} self._forbidden_elements: list[int] | None = None - self._scale_vdw_radii: float = 4.0 / 3.0 - self._scale_minimal_bondlength: float = 0.75 + self._scale_vdw_radii: float = 1.25 + self._scale_minimal_bondlength: float = 0.8 + self._contract_coords: bool = False def get_identifier(self) -> str: return "generate" @@ -383,6 +384,22 @@ def scale_minimal_bondlength(self, scale_minimal_bondlength: float): raise ValueError("Scale minimal bond length should be greater than 0.") self._scale_minimal_bondlength = scale_minimal_bondlength + @property + def contract_coords(self): + """ + Get the contract_coords flag. + """ + return self._contract_coords + + @contract_coords.setter + def contract_coords(self, contract_coords: bool): + """ + Set the contract_coords flag. + """ + if not isinstance(contract_coords, bool): + raise TypeError("Contract coords should be a boolean.") + self._contract_coords = contract_coords + class RefineConfig(BaseConfig): """