Typing: Parameters, vector_algebra

- Parameters typing - vector_algebra tests and refactoring - Various pyright typing improvements
2023-12-12 20:41:17 +01:00
parent 4106726483
commit 57ad5d8384
18 changed files with 488 additions and 183 deletions
--- a/propka/conformation_container.py
+++ b/propka/conformation_container.py
@@ -7,6 +7,8 @@ Container data structure for molecular conformations.
 import logging
 import functools
 from typing import Iterable, List, NoReturn, Optional, TYPE_CHECKING, Set
+from propka.lib import Options
+from propka.version import Version

 if TYPE_CHECKING:
    from propka.atom import Atom
@@ -19,6 +21,7 @@ from propka.coupled_groups import NCCG
 from propka.determinants import set_backbone_determinants, set_ion_determinants
 from propka.determinants import set_determinants
 from propka.group import Group, is_group
+from propka.parameters import Parameters


 _LOGGER = logging.getLogger(__name__)
@@ -44,9 +47,9 @@ class ConformationContainer:
    """

    def __init__(self,
-                 name: str = '',
-                 parameters=None,
-                 molecular_container: Optional["MolecularContainer"] = None):
+                 name: str,
+                 parameters: Parameters,
+                 molecular_container: "MolecularContainer"):
        """Initialize conformation container.

        Args:
@@ -145,6 +148,7 @@ class ConformationContainer:
        Returns:
            a set of bonded atom groups
        """
+        assert self.parameters is not None
        res: Set[Group] = set()
        for bond_atom in atom.bonded_atoms:
            # skip the original atom
@@ -187,7 +191,7 @@ class ConformationContainer:

        # If --titrate_only option is set, make non-specified residues
        # un-titratable:
-        assert self.molecular_container is not None
+        assert self.molecular_container.options is not None
        titrate_only = self.molecular_container.options.titrate_only
        if titrate_only is not None:
            atom = group.atom
@@ -196,7 +200,7 @@ class ConformationContainer:
                if group.residue_type == 'CYS':
                    group.exclude_cys_from_results = True

-    def calculate_pka(self, version, options):
+    def calculate_pka(self, version: Version, options: Options):
        """Calculate pKas for conformation container.

        Args:
@@ -349,7 +353,7 @@ class ConformationContainer:
                                                  reference=reference)
        return ddg

-    def calculate_charge(self, parameters, ph=None):
+    def calculate_charge(self, parameters, ph: float):
        """Calculate charge for folded and unfolded states.

        Args:
--- a/propka/coupled_groups.py
+++ b/propka/coupled_groups.py
@@ -33,6 +33,7 @@ class NonCovalentlyCoupledGroups:
        Returns:
            dictionary describing coupling
        """
+        assert self.parameters is not None
        # check if the interaction energy is high enough
        interaction_energy = max(self.get_interaction(group1, group2),
                                 self.get_interaction(group2, group1))
@@ -105,6 +106,7 @@ class NonCovalentlyCoupledGroups:
        Returns:
            float value of scaling factor
        """
+        assert self.parameters is not None
        intrinsic_pka_diff = abs(pka1-pka2)
        res = 0.0
        if intrinsic_pka_diff <= self.parameters.max_intrinsic_pka_diff:
@@ -122,6 +124,7 @@ class NonCovalentlyCoupledGroups:
        Returns:
            float value of scaling factor
        """
+        assert self.parameters is not None
        free_energy_diff = abs(energy1-energy2)
        res = 0.0
        if free_energy_diff <= self.parameters.max_free_energy_diff:
@@ -136,6 +139,7 @@ class NonCovalentlyCoupledGroups:
        Returns:
            float value of scaling factor
        """
+        assert self.parameters is not None
        res = 0.0
        interaction_energy = abs(interaction_energy)
        if interaction_energy >= self.parameters.min_interaction_energy:
--- a/propka/determinants.py
+++ b/propka/determinants.py
@@ -14,12 +14,18 @@ Functions to manipulate :class:`propka.determinant.Determinant` objects.

 """
 import math
+from typing import List
+
+import propka.calculations
 import propka.iterative
 import propka.lib
 import propka.vector_algebra
 from propka.calculations import squared_distance, get_smallest_distance
 from propka.energy import angle_distance_factors, hydrogen_bond_energy
 from propka.determinant import Determinant
+from propka.group import Group
+from propka.iterative import Interaction
+from propka.version import Version


 # Cutoff for angle factor
@@ -28,7 +34,7 @@ from propka.determinant import Determinant
 FANGLE_MIN = 0.001


-def set_determinants(propka_groups, version=None, options=None):
+def set_determinants(propka_groups: List[Group], version: Version, options=None):
    """Add side-chain and coulomb determinants/perturbations to all residues.

    NOTE - backbone determinants are set separately
@@ -38,7 +44,7 @@ def set_determinants(propka_groups, version=None, options=None):
        version:  version object
        options:  options object
    """
-    iterative_interactions = []
+    iterative_interactions: List[Interaction] = []
    # --- NonIterative section ---#
    for group1 in propka_groups:
        for group2 in propka_groups:
@@ -77,7 +83,7 @@ def add_determinants(group1, group2, distance, version):
    add_coulomb_determinants(group1, group2, distance, version)


-def add_sidechain_determinants(group1, group2, version=None):
+def add_sidechain_determinants(group1: Group, group2: Group, version: Version):
    """Add side-chain determinants and perturbations.

    NOTE - res_num1 > res_num2
@@ -236,6 +242,8 @@ def set_backbone_determinants(titratable_groups, backbone_groups, version):
                get_smallest_distance(
                    backbone_interaction_atoms,
                    titratable_group_interaction_atoms))
+            assert backbone_atom is not None
+            assert titratable_atom is not None
            # get the parameters
            parameters = (
                version.get_backbone_hydrogen_bond_parameters(
--- a/propka/energy.py
+++ b/propka/energy.py
@@ -7,11 +7,15 @@ Energy calculations.
 """
 import math
 import logging
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Optional, Sequence
+
+from propka.atom import Atom
+from propka.parameters import Parameters

 if TYPE_CHECKING:
    from propka.conformation_container import ConformationContainer
    from propka.group import Group
+    from propka.version import Version

 from propka.calculations import squared_distance, get_smallest_distance

@@ -89,7 +93,7 @@ def calculate_scale_factor(parameters, weight: float) -> float:
    return scale_factor


-def calculate_weight(parameters, num_volume: int) -> float:
+def calculate_weight(parameters: Parameters, num_volume: float) -> float:
    """Calculate the atom-based desolvation weight.

    TODO - figure out why a similar function exists in version.py
@@ -109,7 +113,7 @@ def calculate_weight(parameters, num_volume: int) -> float:
    return weight


-def calculate_pair_weight(parameters, num_volume1: int, num_volume2: int) -> float:
+def calculate_pair_weight(parameters: Parameters, num_volume1: int, num_volume2: int) -> float:
    """Calculate the atom-pair based desolvation weight.

    Args:
@@ -148,7 +152,11 @@ def hydrogen_bond_energy(dist, dpka_max: float, cutoffs, f_angle=1.0) -> float:
    return abs(dpka)


-def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None):
+def angle_distance_factors(
+        atom1: Optional[Atom] = None,
+        atom2: Atom = None,  # type: ignore[assignment]
+        atom3: Atom = None,  # type: ignore[assignment]
+        center: Optional[Sequence[float]] = None):
    """Calculate distance and angle factors for three atoms for backbone
    interactions.

@@ -182,6 +190,7 @@ def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None):
    dy_32 = dy_32/dist_23
    dz_32 = dz_32/dist_23
    if atom1 is None:
+        assert center is not None
        dx_21 = center[0] - atom2.x
        dy_21 = center[1] - atom2.y
        dz_21 = center[2] - atom2.z
@@ -197,7 +206,7 @@ def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None):
    return dist_12, f_angle, dist_23


-def hydrogen_bond_interaction(group1, group2, version):
+def hydrogen_bond_interaction(group1: "Group", group2: "Group", version: "Version"):
    """Calculate energy for hydrogen bond interactions between two groups.

    Args:
@@ -213,7 +222,7 @@ def hydrogen_bond_interaction(group1, group2, version):
    [closest_atom1, dist, closest_atom2] = get_smallest_distance(
        atoms1, atoms2
    )
-    if None in [closest_atom1, closest_atom2]:
+    if closest_atom1 is None or closest_atom2 is None:
        _LOGGER.warning(
            'Side chain interaction failed for {0:s} and {1:s}'.format(
                group1.label, group2.label))
@@ -297,7 +306,7 @@ def electrostatic_interaction(group1, group2, dist, version):
    return value


-def check_coulomb_pair(parameters, group1, group2, dist):
+def check_coulomb_pair(parameters: Parameters, group1: "Group", group2: "Group", dist: float) -> bool:
    """Checks if this Coulomb interaction should be done.

    NOTE - this is a propka2.0 hack
--- a/propka/group.py
+++ b/propka/group.py
@@ -14,6 +14,7 @@ import math
 from typing import cast, Dict, Iterable, List, NoReturn, Optional

 import propka.ligand
+from propka.parameters import Parameters
 import propka.protonate
 from propka.atom import Atom
 from propka.ligand_pka_values import LigandPkaValues
@@ -90,7 +91,7 @@ class Group:
        self.y = 0.0
        self.z = 0.0
        self.charge = 0
-        self.parameters = None
+        self.parameters: Optional[Parameters] = None
        self.exclude_cys_from_results = False
        self.interaction_atoms_for_acids: List[Atom] = []
        self.interaction_atoms_for_bases: List[Atom] = []
@@ -320,6 +321,7 @@ class Group:

    def setup(self):
        """Set up a group."""
+        assert self.parameters is not None
        # set the charges
        if self.type in self.parameters.charge.keys():
            self.charge = self.parameters.charge[self.type]
@@ -403,7 +405,7 @@ class Group:
                    '             {0:s}'.format(
                        str(self.interaction_atoms_for_bases[i])))

-    def get_interaction_atoms(self, interacting_group) -> List[Atom]:
+    def get_interaction_atoms(self, interacting_group: "Group") -> List[Atom]:
        """Get atoms involved in interaction with other group.

        Args:
@@ -592,7 +594,7 @@ class Group:
        ddg = ddg_neutral + ddg_low
        return ddg

-    def calculate_charge(self, _, ph=7.0, state='folded'):
+    def calculate_charge(self, _, ph: float = 7.0, state: str = 'folded') -> float:
        """Calculate the charge of the specified state at the specified pH.

        Args:
@@ -610,7 +612,7 @@ class Group:
        charge = self.charge*(conc_ratio/(1.0+conc_ratio))
        return charge

-    def use_in_calculations(self):
+    def use_in_calculations(self) -> bool:
        """Indicate whether group should be included in results report.

        If --titrate_only option is specified, only residues that are
@@ -1219,7 +1221,7 @@ class TitratableLigandGroup(Group):
        self.model_pka_set = True


-def is_group(parameters, atom: Atom) -> Optional[Group]:
+def is_group(parameters: Parameters, atom: Atom) -> Optional[Group]:
    """Identify whether the atom belongs to a group.

    Args:
@@ -1246,7 +1248,7 @@ def is_group(parameters, atom: Atom) -> Optional[Group]:
        ligand_group = is_ligand_group_by_groups(parameters, atom)
    else:
        raise Exception(
-            'Unknown ligand typing method \'{0.s}\''.format(
+            'Unknown ligand typing method \'{0:s}\''.format(
                parameters.ligand_typing))
    if ligand_group:
        return ligand_group
@@ -1369,7 +1371,7 @@ def is_ligand_group_by_groups(_, atom: Atom) -> Optional[Group]:
    return None


-def is_ligand_group_by_marvin_pkas(parameters, atom: Atom) -> Optional[Group]:
+def is_ligand_group_by_marvin_pkas(parameters: Parameters, atom: Atom) -> Optional[Group]:
    """Identify whether the atom belongs to a ligand group by calculating
    'Marvin pKas'.

--- a/propka/hydrogens.py
+++ b/propka/hydrogens.py
@@ -122,6 +122,12 @@ def add_arg_hydrogen(residue: List[Atom]) -> List[Atom]:
    Returns:
        list of hydrogen atoms
    """
+    cd_atom: Optional[Atom] = None
+    cz_atom: Optional[Atom] = None
+    ne_atom: Optional[Atom] = None
+    nh1_atom: Optional[Atom] = None
+    nh2_atom: Optional[Atom] = None
+
    for atom in residue:
        if atom.name == "CD":
            cd_atom = atom
@@ -133,6 +139,11 @@ def add_arg_hydrogen(residue: List[Atom]) -> List[Atom]:
            nh1_atom = atom
        elif atom.name == "NH2":
            nh2_atom = atom
+
+    if (cd_atom is None or cz_atom is None or ne_atom is None or nh1_atom is None
+            or nh2_atom is None):
+        raise ValueError("Unable to find all atoms")
+
    h1_atom = protonate_sp2(cd_atom, ne_atom, cz_atom)
    h1_atom.name = "HE"
    h2_atom = protonate_direction(nh1_atom, ne_atom, cz_atom)
@@ -152,6 +163,12 @@ def add_his_hydrogen(residue: List[Atom]) -> None:
    Args:
        residue:  histidine residue to protonate
    """
+    cg_atom: Optional[Atom] = None
+    nd_atom: Optional[Atom] = None
+    cd_atom: Optional[Atom] = None
+    ce_atom: Optional[Atom] = None
+    ne_atom: Optional[Atom] = None
+
    for atom in residue:
        if atom.name == "CG":
            cg_atom = atom
@@ -163,6 +180,11 @@ def add_his_hydrogen(residue: List[Atom]) -> None:
            ce_atom = atom
        elif atom.name == "NE2":
            ne_atom = atom
+
+    if (cg_atom is None or nd_atom is None or cd_atom is None or ce_atom is None
+            or ne_atom is None):
+        raise ValueError("Unable to find all atoms")
+
    hd_atom = protonate_sp2(cg_atom, nd_atom, ce_atom)
    hd_atom.name = "HND"
    he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom)
@@ -177,6 +199,7 @@ def add_trp_hydrogen(residue: List[Atom]) -> None:
    """
    cd_atom = None
    ne_atom = None
+    ce_atom = None
    for atom in residue:
        if atom.name == "CD1":
            cd_atom = atom
--- a/propka/iterative.py
+++ b/propka/iterative.py
@@ -7,7 +7,10 @@ involve :class:`propka.determinant.Determinant` instances.

 """
 import logging
+from typing import Dict, Iterable, List, Optional, Sequence, Tuple
 from propka.determinant import Determinant
+from propka.group import Group
+from propka.version import Version


 _LOGGER = logging.getLogger(__name__)
@@ -16,9 +19,11 @@ _LOGGER = logging.getLogger(__name__)
 # TODO - these are undocumented constants
 UNK_MIN_VALUE = 0.005

+Interaction = list

-def add_to_determinant_list(group1, group2, distance, iterative_interactions,
-                            version):
+
+def add_to_determinant_list(group1: Group, group2: Group, distance: float,
+                            iterative_interactions: List[Interaction], version: Version):
    """Add iterative determinantes to the list.

    [[R1, R2], [side-chain, coulomb], [A1, A2]], ...
@@ -48,7 +53,8 @@ def add_to_determinant_list(group1, group2, distance, iterative_interactions,
        iterative_interactions.append(interaction)


-def add_iterative_acid_pair(object1, object2, interaction):
+def add_iterative_acid_pair(object1: "Iterative", object2: "Iterative",
+                            interaction: Interaction):
    """Add the Coulomb 'iterative' interaction (an acid pair).

    The higher pKa is raised  with QQ+HB
@@ -90,7 +96,8 @@ def add_iterative_acid_pair(object1, object2, interaction):
        annihilation[1] = -diff


-def add_iterative_base_pair(object1, object2, interaction):
+def add_iterative_base_pair(object1: "Iterative", object2: "Iterative",
+                            interaction: Interaction):
    """Add the Coulomb 'iterative' interaction (a base pair).

    The lower pKa is lowered
@@ -132,7 +139,8 @@ def add_iterative_base_pair(object1, object2, interaction):
        annihilation[1] = -diff


-def add_iterative_ion_pair(object1, object2, interaction, version):
+def add_iterative_ion_pair(object1: "Iterative", object2: "Iterative",
+                           interaction: Interaction, version: Version):
    """Add the Coulomb 'iterative' interaction (an acid-base pair)

    the pKa of the acid is lowered & the pKa of the base is raised
@@ -194,7 +202,7 @@ def add_iterative_ion_pair(object1, object2, interaction, version):
                object2.determinants['sidechain'].append(interaction)


-def add_determinants(iterative_interactions, version, _=None):
+def add_determinants(iterative_interactions: List[Interaction], version: Version, _=None):
    """Add determinants iteratively.

    The iterative pKa scheme. Later it is all added in 'calculateTotalPKA'
@@ -205,7 +213,7 @@ def add_determinants(iterative_interactions, version, _=None):
        _:  options object
    """
    # --- setup ---
-    iteratives = []
+    iteratives: List[Iterative] = []
    done_group = []
    # create iterative objects with references to their real group counterparts
    for interaction in iterative_interactions:
@@ -270,6 +278,7 @@ def add_determinants(iterative_interactions, version, _=None):

        # reset pka_old & storing pka_new in pka_iter
        for itres in iteratives:
+            assert itres.pka_new is not None
            itres.pka_old = itres.pka_new
            itres.pka_iter.append(itres.pka_new)

@@ -302,7 +311,10 @@ def add_determinants(iterative_interactions, version, _=None):
                    itres.group.determinants[type_].append(new_det)


-def find_iterative(pair, iteratives):
+def find_iterative(
+    pair: Sequence[Group],
+    iteratives: Iterable["Iterative"],
+) -> Tuple["Iterative", "Iterative"]:
    """Find the 'iteratives' that correspond to the groups in 'pair'.

    Args:
@@ -312,11 +324,15 @@ def find_iterative(pair, iteratives):
        1. first matched iterative
        2. second matched iterative
    """
+    iterative0: Optional[Iterative] = None
+    iterative1: Optional[Iterative] = None
    for iterative in iteratives:
        if iterative.group == pair[0]:
            iterative0 = iterative
        elif iterative.group == pair[1]:
            iterative1 = iterative
+    if iterative0 is None or iterative1 is None:
+        raise LookupError("iteratives not found")
    return iterative0, iterative1


@@ -327,7 +343,7 @@ class Iterative:
    after the iterations are finished.
    """

-    def __init__(self, group):
+    def __init__(self, group: Group):
        """Initialize object with group.

        Args:
@@ -337,11 +353,15 @@ class Iterative:
        self.atom = group.atom
        self.res_name = group.residue_type
        self.q = group.charge
-        self.pka_old = None
-        self.pka_new = None
-        self.pka_iter = []
+        self.pka_old: Optional[float] = None
+        self.pka_new: Optional[float] = None
+        self.pka_iter: List[float] = []
        self.pka_noniterative = 0.00
-        self.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []}
+        self.determinants: Dict[str, list] = {
+            'sidechain': [],
+            'backbone': [],
+            'coulomb': []
+        }
        self.group = group
        self.converged = True
        # Calculate the Non-Iterative part of pKa from the group object
@@ -368,8 +388,9 @@ class Iterative:
        self.pka_noniterative += coulomb
        self.pka_old = self.pka_noniterative

-    def __eq__(self, other):
+    def __eq__(self, other) -> bool:
        """Needed to use objects in sets."""
+        assert isinstance(other, (Iterative, Group)), type(other)
        if self.atom.type == 'atom':
            # In case of protein atoms we trust the labels
            return self.label == other.label
--- a/propka/lib.py
+++ b/propka/lib.py
@@ -9,6 +9,10 @@ import sys
 import logging
 import argparse
 from pathlib import Path
+from typing import List, TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from propka.atom import Atom


 _LOGGER = logging.getLogger(__name__)
@@ -19,6 +23,8 @@ EXPECTED_ATOM_NUMBERS = {'ALA': 5, 'ARG': 11, 'ASN': 8, 'ASP': 8, 'CYS': 6,
                         'LEU': 8, 'LYS': 9, 'MET': 8, 'PHE': 11, 'PRO': 7,
                         'SER': 6, 'THR': 7, 'TRP': 14, 'TYR': 12, 'VAL': 7}

+Options = argparse.Namespace
+

 def protein_precheck(conformations, names):
    """Check protein for correct number of atoms, etc.
@@ -73,7 +79,7 @@ def resid_from_atom(atom):
        atom.res_num, atom.chain_id, atom.icode)


-def split_atoms_into_molecules(atoms):
+def split_atoms_into_molecules(atoms: List["Atom"]):
    """Maps atoms into molecules.

    Args:
@@ -81,14 +87,14 @@ def split_atoms_into_molecules(atoms):
    Returns:
        list of molecules
    """
-    molecules = []
+    molecules: List[List["Atom"]] = []
    while len(atoms) > 0:
        initial_atom = atoms.pop()
        molecules.append(make_molecule(initial_atom, atoms))
    return molecules


-def make_molecule(atom, atoms):
+def make_molecule(atom: "Atom", atoms: List["Atom"]):
    """Make a molecule from atoms.

    Args:
@@ -302,7 +308,7 @@ def build_parser(parser=None):
    return parser


-def loadOptions(args=None):
+def loadOptions(args=None) -> Options:
    """
    Load the arguments parser with options. Note that verbosity is set as soon
    as this function is invoked.
@@ -315,7 +321,7 @@ def loadOptions(args=None):
    # loading the parser
    parser = build_parser()
    # parsing and returning options and arguments
-    options = parser.parse_args(args)
+    options = parser.parse_args(args, namespace=Options())

    # adding specified filenames to arguments
    options.filenames.append(options.input_pdb)
--- a/propka/ligand_pka_values.py
+++ b/propka/ligand_pka_values.py
@@ -11,11 +11,13 @@ programs are required).
 """
 import logging
 import os
+import shutil
 import subprocess
 import sys
 import warnings
 from propka.output import write_mol2_for_atoms
 from propka.lib import split_atoms_into_molecules
+from propka.parameters import Parameters


 _LOGGER = logging.getLogger(__name__)
@@ -24,7 +26,7 @@ _LOGGER = logging.getLogger(__name__)
 class LigandPkaValues:
    """Ligand pKa value class."""

-    def __init__(self, parameters):
+    def __init__(self, parameters: Parameters):
        """Initialize object with parameters.

        Args:
@@ -47,20 +49,16 @@ class LigandPkaValues:
        Returns:
            location of program
        """
-        path = os.environ.get('PATH').split(os.pathsep)
-        locs = [
-            i for i in filter(lambda loc: os.access(loc, os.F_OK),
-                              map(lambda dir: os.path.join(dir, program),
-                                  path))]
-        if len(locs) == 0:
+        loc = shutil.which(program)
+        if loc is None:
            str_ = "'Error: Could not find {0:s}.".format(program)
            str_ += ' Please make sure that it is found in the path.'
            _LOGGER.info(str_)
            sys.exit(-1)
-        return locs[0]
+        return loc

    def get_marvin_pkas_for_pdb_file(
-            self, molecule, parameters, num_pkas=10, min_ph=-10, max_ph=20):
+            self, molecule, parameters, num_pkas=10, min_ph=-10.0, max_ph=20.0):
        """Use Marvin executables to get pKas for a PDB file.

        Args:
@@ -75,7 +73,7 @@ class LigandPkaValues:
            molecule, num_pkas=num_pkas, min_ph=min_ph, max_ph=max_ph)

    def get_marvin_pkas_for_molecular_container(self, molecule, num_pkas=10,
-                                                min_ph=-10, max_ph=20):
+                                                min_ph=-10.0, max_ph=20.0):
        """Use Marvin executables to calculate pKas for a molecular container.

        Args:
@@ -93,8 +91,8 @@ class LigandPkaValues:

    def get_marvin_pkas_for_conformation_container(self, conformation,
                                                   name='temp', reuse=False,
-                                                   num_pkas=10, min_ph=-10,
-                                                   max_ph=20):
+                                                   num_pkas=10, min_ph=-10.0,
+                                                   max_ph=20.0):
        """Use Marvin executables to calculate pKas for a conformation container.

        Args:
@@ -111,7 +109,7 @@ class LigandPkaValues:
            num_pkas=num_pkas, min_ph=min_ph, max_ph=max_ph)

    def get_marvin_pkas_for_atoms(self, atoms, name='temp', reuse=False,
-                                  num_pkas=10, min_ph=-10, max_ph=20):
+                                  num_pkas=10, min_ph=-10.0, max_ph=20.0):
        """Use Marvin executables to calculate pKas for a list of atoms.

        Args:
@@ -131,8 +129,8 @@ class LigandPkaValues:
                min_ph=min_ph, max_ph=max_ph)

    def get_marvin_pkas_for_molecule(self, atoms, filename='__tmp_ligand.mol2',
-                                     reuse=False, num_pkas=10, min_ph=-10,
-                                     max_ph=20):
+                                     reuse=False, num_pkas=10, min_ph=-10.0,
+                                     max_ph=20.0):
        """Use Marvin executables to calculate pKas for a molecule.

        Args:
--- a/propka/molecular_container.py
+++ b/propka/molecular_container.py
@@ -7,11 +7,12 @@ Molecular container for storing all contents of PDB files.
 import logging
 import os
 from typing import Dict, List, Optional, Tuple
+from propka.parameters import Parameters

 import propka.version
 from propka.output import write_pka, print_header, print_result
 from propka.conformation_container import ConformationContainer
-from propka.lib import make_grid
+from propka.lib import make_grid, Options


 _LOGGER = logging.getLogger(__name__)
@@ -35,7 +36,7 @@ class MolecularContainer:
    name: Optional[str]
    version: propka.version.Version

-    def __init__(self, parameters, options=None) -> None:
+    def __init__(self, parameters: Parameters, options: Options) -> None:
        """Initialize molecular container.

        Args:
--- a/propka/output.py
+++ b/propka/output.py
@@ -365,7 +365,7 @@ def write_jackal_scap_file(mutation_data=None, filename="1xxx_scap.list",
        for chain_id, _, res_num, code2 in mutation_data:
            str_ = "{chain:s}, {num:d}, {code:s}\n".format(
                chain=chain_id, num=res_num, code=code2)
-        file_.write(str_)
+            file_.write(str_)


 def write_scwrl_sequence_file(sequence, filename="x-ray.seq", _=None):
--- a/propka/parameters.py
+++ b/propka/parameters.py
@@ -11,67 +11,128 @@ in configuration file.

 """
 import logging
+from dataclasses import dataclass, field
+from typing import Dict, List
+
+try:
+    # New in version 3.10, deprecated since version 3.12
+    from typing import TypeAlias
+except ImportError:
+    TypeAlias = "TypeAlias"  # type: ignore


 _LOGGER = logging.getLogger(__name__)


-#: matrices
-MATRICES = ['interaction_matrix']
-#: pari-wise matrices
-PAIR_WISE_MATRICES = ['sidechain_cutoffs']
-#: :class:`dict` containing numbers
-NUMBER_DICTIONARIES = [
-    'VanDerWaalsVolume', 'charge', 'model_pkas', 'ions', 'valence_electrons',
-    'custom_model_pkas']
-#: :class:`dict` containing lists
-LIST_DICTIONARIES = ['backbone_NH_hydrogen_bond', 'backbone_CO_hydrogen_bond']
-#: :class:`dict` containing strings
-STRING_DICTIONARIES = ['protein_group_mapping']
-#: :class:`list` containing strings
-STRING_LISTS = [
-    'ignore_residues', 'angular_dependent_sidechain_interactions',
-    'acid_list', 'base_list', 'exclude_sidechain_interactions',
-    'backbone_reorganisation_list', 'write_out_order']
-#: distances (:class:`float`)
-DISTANCES = ['desolv_cutoff', 'buried_cutoff', 'coulomb_cutoff1',
-             'coulomb_cutoff2']
-#: other parameters
-PARAMETERS = [
-    'Nmin', 'Nmax', 'desolvationSurfaceScalingFactor', 'desolvationPrefactor',
-    'desolvationAllowance', 'coulomb_diel', 'COO_HIS_exception',
-    'OCO_HIS_exception', 'CYS_HIS_exception', 'CYS_CYS_exception',
-    'min_ligand_model_pka', 'max_ligand_model_pka',
-    'include_H_in_interactions', 'coupling_max_number_of_bonds',
-    'min_bond_distance_for_hydrogen_bonds', 'coupling_penalty',
-    'shared_determinants', 'common_charge_centre', 'hide_penalised_group',
-    'remove_penalised_group', 'max_intrinsic_pka_diff',
-    'min_interaction_energy', 'max_free_energy_diff', 'min_swap_pka_shift',
-    'min_pka', 'max_pka', 'sidechain_interaction']
-# :class:`str` parameters
-STRINGS = ['version', 'output_file_tag', 'ligand_typing', 'pH', 'reference']
+class squared_property:
+
+    def __set_name__(self, owner, name: str):
+        assert name.endswith("_squared")
+        self._name_not_squared = name[:-len("_squared")]  # removesuffix()
+
+    def __get__(self, instance, owner=None) -> float:
+        if instance is None:
+            return self  # type: ignore[return-value]
+        return getattr(instance, self._name_not_squared)**2
+
+    def __set__(self, instance, value: float):
+        setattr(instance, self._name_not_squared, value**0.5)


+_T_MATRIX: TypeAlias = "InteractionMatrix"
+_T_PAIR_WISE_MATRIX: TypeAlias = "PairwiseMatrix"
+_T_NUMBER_DICTIONARY = Dict[str, float]
+_T_LIST_DICTIONARY = Dict[str, list]
+_T_STRING_DICTIONARY = Dict[str, str]
+_T_STRING_LIST = List[str]
+_T_STRING = str
+_T_BOOL = bool
+
+
+@dataclass
 class Parameters:
    """PROPKA parameter class."""

-    def __init__(self):
-        """Initialize parameter class.
+    # MATRICES
+    interaction_matrix: _T_MATRIX = field(
+        default_factory=lambda: InteractionMatrix("interaction_matrix"))

-        Args:
-            parameter_file:  file with parameters
-        """
-        # TODO - need to define all members explicitly
-        self.model_pkas = {}
-        self.interaction_matrix = InteractionMatrix("interaction_matrix")
-        self.sidechain_cutoffs = None
-        # TODO - it would be nice to rename these; they're defined everywhere
-        self.COO_HIS_exception = None
-        self.OCO_HIS_exception = None
-        self.CYS_HIS_exception = None
-        self.CYS_CYS_exception = None
-        # These functions set up remaining data structures implicitly
-        self.set_up_data_structures()
+    # PAIR_WISE_MATRICES
+    sidechain_cutoffs: _T_PAIR_WISE_MATRIX = field(
+        default_factory=lambda: PairwiseMatrix("sidechain_cutoffs"))
+
+    # NUMBER_DICTIONARIES
+    VanDerWaalsVolume: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+    charge: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+    model_pkas: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+    ions: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+    valence_electrons: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+    custom_model_pkas: _T_NUMBER_DICTIONARY = field(default_factory=dict)
+
+    # LIST_DICTIONARIES
+    backbone_NH_hydrogen_bond: _T_LIST_DICTIONARY = field(default_factory=dict)
+    backbone_CO_hydrogen_bond: _T_LIST_DICTIONARY = field(default_factory=dict)
+
+    # STRING_DICTIONARIES
+    protein_group_mapping: _T_STRING_DICTIONARY = field(default_factory=dict)
+
+    # STRING_LISTS
+    ignore_residues: _T_STRING_LIST = field(default_factory=list)
+    angular_dependent_sidechain_interactions: _T_STRING_LIST = field(default_factory=list)
+    acid_list: _T_STRING_LIST = field(default_factory=list)
+    base_list: _T_STRING_LIST = field(default_factory=list)
+    exclude_sidechain_interactions: _T_STRING_LIST = field(default_factory=list)
+    backbone_reorganisation_list: _T_STRING_LIST = field(default_factory=list)
+    write_out_order: _T_STRING_LIST = field(default_factory=list)
+
+    # DISTANCES
+    desolv_cutoff: float = 20.0
+    buried_cutoff: float = 15.0
+    coulomb_cutoff1: float = 4.0
+    coulomb_cutoff2: float = 10.0
+
+    # DISTANCES SQUARED
+    desolv_cutoff_squared = squared_property()
+    buried_cutoff_squared = squared_property()
+    coulomb_cutoff1_squared = squared_property()
+    coulomb_cutoff2_squared = squared_property()
+
+    # STRINGS
+    version: _T_STRING = "VersionA"
+    output_file_tag: _T_STRING = ""
+    ligand_typing: _T_STRING = "groups"
+    pH: _T_STRING = "variable"
+    reference: _T_STRING = "neutral"
+
+    # PARAMETERS
+    Nmin: int = 280
+    Nmax: int = 560
+    desolvationSurfaceScalingFactor: float = 0.25
+    desolvationPrefactor: float = -13.0
+    desolvationAllowance: float = 0.0
+    coulomb_diel: float = 80.0
+    # TODO - it would be nice to rename these; they're defined everywhere
+    COO_HIS_exception: float = 1.60
+    OCO_HIS_exception: float = 1.60
+    CYS_HIS_exception: float = 1.60
+    CYS_CYS_exception: float = 3.60
+    min_ligand_model_pka: float = -10.0
+    max_ligand_model_pka: float = 20.0
+    # include_H_in_interactions: NoReturn = None
+    coupling_max_number_of_bonds: int = 3
+    min_bond_distance_for_hydrogen_bonds: int = 4
+    # coupling_penalty: NoReturn = None
+    shared_determinants: _T_BOOL = False
+    common_charge_centre: _T_BOOL = False
+    # hide_penalised_group: NoReturn = None
+    remove_penalised_group: _T_BOOL = True
+    max_intrinsic_pka_diff: float = 2.0
+    min_interaction_energy: float = 0.5
+    max_free_energy_diff: float = 1.0
+    min_swap_pka_shift: float = 1.0
+    min_pka: float = 0.0
+    max_pka: float = 10.0
+    sidechain_interaction: float = 0.85

    def parse_line(self, line):
        """Parse parameter file line."""
@@ -84,22 +145,21 @@ class Parameters:
        if len(words) == 0:
            return
        # parse the words
-        if len(words) == 3 and words[0] in NUMBER_DICTIONARIES:
+        typeannotation = self.__annotations__.get(words[0])
+        if typeannotation is _T_NUMBER_DICTIONARY:
            self.parse_to_number_dictionary(words)
-        elif len(words) == 2 and words[0] in STRING_LISTS:
+        elif typeannotation is _T_STRING_LIST:
            self.parse_to_string_list(words)
-        elif len(words) == 2 and words[0] in DISTANCES:
-            self.parse_distance(words)
-        elif len(words) == 2 and words[0] in PARAMETERS:
-            self.parse_parameter(words)
-        elif len(words) == 2 and words[0] in STRINGS:
+        elif typeannotation is _T_STRING:
            self.parse_string(words)
-        elif len(words) > 2 and words[0] in LIST_DICTIONARIES:
+        elif typeannotation is _T_LIST_DICTIONARY:
            self.parse_to_list_dictionary(words)
-        elif words[0] in MATRICES+PAIR_WISE_MATRICES:
+        elif typeannotation is _T_MATRIX or typeannotation is _T_PAIR_WISE_MATRIX:
            self.parse_to_matrix(words)
-        elif len(words) == 3 and words[0] in STRING_DICTIONARIES:
+        elif typeannotation is _T_STRING_DICTIONARY:
            self.parse_to_string_dictionary(words)
+        else:
+            self.parse_parameter(words)

    def parse_to_number_dictionary(self, words):
        """Parse field to number dictionary.
@@ -107,6 +167,7 @@ class Parameters:
        Args:
            words:  strings to parse.
        """
+        assert len(words) == 3, words
        dict_ = getattr(self, words[0])
        key = words[1]
        value = words[2]
@@ -118,17 +179,19 @@ class Parameters:
        Args:
            words:  strings to parse
        """
+        assert len(words) == 3, words
        dict_ = getattr(self, words[0])
        key = words[1]
        value = words[2]
        dict_[key] = value

-    def parse_to_list_dictionary(self, words):
+    def parse_to_list_dictionary(self, words: List[str]):
        """Parse field to list dictionary.

        Args:
            words:  strings to parse.
        """
+        assert len(words) > 2, words
        dict_ = getattr(self, words[0])
        key = words[1]
        if key not in dict_:
@@ -144,6 +207,7 @@ class Parameters:
        Args:
            words:  strings to parse
        """
+        assert len(words) == 2, words
        list_ = getattr(self, words[0])
        value = words[1]
        list_.append(value)
@@ -158,24 +222,13 @@ class Parameters:
        value = tuple(words[1:])
        matrix.add(value)

-    def parse_distance(self, words):
-        """Parse field to distance.
-
-        Args:
-            words:  strings to parse
-        """
-        value = float(words[1])
-        setattr(self, words[0], value)
-        value_sq = value*value
-        attr = "{0:s}_squared".format(words[0])
-        setattr(self, attr, value_sq)
-
    def parse_parameter(self, words):
        """Parse field to parameters.

        Args:
            words:  strings to parse
        """
+        assert len(words) == 2, words
        value = float(words[1])
        setattr(self, words[0], value)

@@ -185,28 +238,9 @@ class Parameters:
        Args:
            words:  strings to parse
        """
+        assert len(words) == 2, words
        setattr(self, words[0], words[1])

-    def set_up_data_structures(self):
-        """Set up internal data structures.
-
-        TODO - it would be better to make these assignments explicit in
-        __init__.
-        """
-        for key_word in (NUMBER_DICTIONARIES + LIST_DICTIONARIES
-                         + STRING_DICTIONARIES):
-            setattr(self, key_word, {})
-        for key_word in STRING_LISTS:
-            setattr(self, key_word, [])
-        for key_word in STRINGS:
-            setattr(self, key_word, "")
-        for key_word in MATRICES:
-            matrix = InteractionMatrix(key_word)
-            setattr(self, key_word, matrix)
-        for key_word in PAIR_WISE_MATRICES:
-            matrix = PairwiseMatrix(key_word)
-            setattr(self, key_word, matrix)
-
    def print_interaction_parameters(self):
        """Print interaction parameters."""
        _LOGGER.info('--------------- Model pKa values ----------------------')
--- a/propka/propka.cfg
+++ b/propka/propka.cfg
@@ -107,7 +107,7 @@ interaction_matrix SH  I N N N N N N N N N N I I I I I I N N N N N N N N N N N I
 sidechain_cutoffs default 3.0  4.0  
 # COO
 sidechain_cutoffs COO COO 2.5  3.5
-Sidechain_cutoffs COO SER 2.65 3.65
+sidechain_cutoffs COO SER 2.65 3.65
 sidechain_cutoffs COO ARG 1.85 2.85
 sidechain_cutoffs COO LYS 2.85 3.85
 sidechain_cutoffs COO HIS 2.0  3.0
--- a/propka/protonate.py
+++ b/propka/protonate.py
@@ -9,10 +9,16 @@ protons.
 """
 import logging
 import math
+from typing import Iterable, TYPE_CHECKING
+
 import propka.bonds
 import propka.atom
+from propka.atom import Atom
 from propka.vector_algebra import rotate_vector_around_an_axis, Vector

+if TYPE_CHECKING:
+    from propka.molecular_container import MolecularContainer
+

 _LOGGER = logging.getLogger(__name__)

@@ -58,7 +64,7 @@ class Protonate:
            'Br': 1.41, 'I': 1.61, 'S': 1.35}
        self.protonation_methods = {4: self.tetrahedral, 3: self.trigonal}

-    def protonate(self, molecules):
+    def protonate(self, molecules: "MolecularContainer"):
        """Protonate all atoms in the molecular container.

        Args:
@@ -75,7 +81,7 @@ class Protonate:
                self.protonate_atom(atom)

    @staticmethod
-    def remove_all_hydrogen_atoms(molecular_container):
+    def remove_all_hydrogen_atoms(molecular_container: "MolecularContainer"):
        """Remove all hydrogen atoms from molecule.

        Args:
@@ -86,7 +92,7 @@ class Protonate:
                molecular_container.conformations[name]
                .get_non_hydrogen_atoms())

-    def set_charge(self, atom):
+    def set_charge(self, atom: Atom):
        """Set charge for atom.

        Args:
@@ -109,7 +115,7 @@ class Protonate:
                atom.sybyl_type = atom.sybyl_type.replace('-', '')
                atom.charge_set = True

-    def protonate_atom(self, atom):
+    def protonate_atom(self, atom: Atom):
        """Protonate an atom.

        Args:
@@ -126,7 +132,7 @@ class Protonate:
        atom.is_protonated = True

    @staticmethod
-    def set_proton_names(heavy_atoms):
+    def set_proton_names(heavy_atoms: Iterable[Atom]):
        """Set names for protons.

        Args:
@@ -139,7 +145,7 @@ class Protonate:
                    bonded.name += str(i)
                    i += 1

-    def set_number_of_protons_to_add(self, atom):
+    def set_number_of_protons_to_add(self, atom: Atom):
        """Set the number of protons to add to this atom.

        Args:
@@ -169,7 +175,7 @@ class Protonate:
        _LOGGER.debug('-'*10)
        _LOGGER.debug(atom.number_of_protons_to_add)

-    def set_steric_number_and_lone_pairs(self, atom):
+    def set_steric_number_and_lone_pairs(self, atom: Atom):
        """Set steric number and lone pairs for atom.

        Args:
@@ -207,8 +213,7 @@ class Protonate:
        atom.steric_number += 0
        _LOGGER.debug('{0:>65s}: {1:>4.1f}'.format(
            'Charge(-)', atom.charge))
-        atom.steric_number -= atom.charge
-        atom.steric_number = math.floor(atom.steric_number/2.0)
+        atom.steric_number = math.floor((atom.steric_number - atom.charge) / 2)
        atom.number_of_lone_pairs = (
            atom.steric_number - len(atom.bonded_atoms)
            - atom.number_of_protons_to_add
@@ -220,7 +225,7 @@ class Protonate:
            'Number of lone pairs', atom.number_of_lone_pairs))
        atom.steric_num_lone_pairs_set = True

-    def add_protons(self, atom):
+    def add_protons(self, atom: Atom):
        """Add protons to atom.

        Args:
@@ -236,7 +241,7 @@ class Protonate:
                '(steric number: {0:d})'.format(atom.steric_number)
            )

-    def trigonal(self, atom):
+    def trigonal(self, atom: Atom):
        """Add hydrogens in trigonal geometry.

        Args:
@@ -296,7 +301,7 @@ class Protonate:
            new_a = self.set_bond_distance(new_a, atom.element)
            self.add_proton(atom, cvec+new_a)

-    def tetrahedral(self, atom):
+    def tetrahedral(self, atom: Atom):
        """Protonate atom in tetrahedral geometry.

        Args:
@@ -338,13 +343,14 @@ class Protonate:
            self.add_proton(atom, cvec+new_a)

    @staticmethod
-    def add_proton(atom, position):
+    def add_proton(atom: Atom, position: Vector):
        """Add a proton to an atom at a specific position.

        Args:
            atom:  atom to protonate
            position:  position for proton
        """
+        assert atom.conformation_container is not None
        # Create the new proton
        new_h = propka.atom.Atom()
        new_h.set_property(
@@ -368,7 +374,6 @@ class Protonate:
        new_h.number_of_lone_pairs = 0
        new_h.number_of_protons_to_add = 0
        new_h.num_pi_elec_2_3_bonds = 0
-        new_h.is_protonates = True
        atom.bonded_atoms.append(new_h)
        atom.number_of_protons_to_add -= 1
        atom.conformation_container.add_atom(new_h)
@@ -386,7 +391,7 @@ class Protonate:
                i += 1
        _LOGGER.debug('added %s %s %s', new_h, 'to', atom)

-    def set_bond_distance(self, bvec, element):
+    def set_bond_distance(self, bvec: Vector, element: str) -> Vector:
        """Set bond distance between atom and element.

        Args:
--- a/propka/vector_algebra.py
+++ b/propka/vector_algebra.py
@@ -6,7 +6,7 @@ Vector algebra for PROPKA.
 """
 import logging
 import math
-from typing import Optional, Protocol, Union
+from typing import Optional, Protocol, overload


 _LOGGER = logging.getLogger(__name__)
@@ -69,20 +69,30 @@ class Vector:
                      self.y - other.y,
                      self.z - other.z)

-    def __mul__(self, other: Union["Vector", "Matrix4x4", float]):
+    def dot(self, other: _XYZ) -> float:
+        return self.x * other.x + self.y * other.y + self.z * other.z
+
+    @overload
+    def __mul__(self, other: "Vector") -> float:
+        ...
+
+    @overload
+    def __mul__(self, other: "Matrix4x4") -> "Vector":
+        ...
+
+    @overload
+    def __mul__(self, other: float) -> "Vector":
+        ...
+
+    def __mul__(self, other):
        """Dot product, scalar and matrix multiplication."""
        if isinstance(other, Vector):
-            return self.x * other.x + self.y * other.y + self.z * other.z
-        elif isinstance(other, Matrix4x4):
-            return Vector(
-                xi=other.a11*self.x + other.a12*self.y + other.a13*self.z
-                + other.a14*1.0,
-                yi=other.a21*self.x + other.a22*self.y + other.a23*self.z
-                + other.a24*1.0,
-                zi=other.a31*self.x + other.a32*self.y + other.a33*self.z
-                + other.a34*1.0
-                )
-        elif type(other) in [int, float]:
+            # TODO deprecate in favor of self.dot()
+            return self.dot(other)
+        if isinstance(other, Matrix4x4):
+            # TODO deprecate in favor of matmul operator
+            return other @ self
+        if isinstance(other, (int, float)):
            return Vector(self.x * other, self.y * other, self.z * other)
        raise TypeError(f'{type(other)} not supported')

@@ -90,6 +100,10 @@ class Vector:
        return self.__mul__(other)

    def __pow__(self, other: _XYZ):
+        # TODO deprecate in favor of self.cross()
+        return self.cross(other)
+
+    def cross(self, other: _XYZ):
        """Cross product."""
        return Vector(self.y * other.z - self.z * other.y,
                      self.z * other.x - self.x * other.z,
@@ -160,6 +174,17 @@ class Matrix4x4:
        self.a43 = a43i
        self.a44 = a44i

+    def __matmul__(self, v: _XYZ) -> Vector:
+        """Matrix vector multiplication with homogeneous coordinates.
+
+        Assumes that the last row is (0, 0, 0, 1).
+        """
+        return Vector(
+            self.a11 * v.x + self.a12 * v.y + self.a13 * v.z + self.a14,
+            self.a21 * v.x + self.a22 * v.y + self.a23 * v.z + self.a24,
+            self.a31 * v.x + self.a32 * v.y + self.a33 * v.z + self.a34,
+        )
+

 def angle(avec: Vector, bvec: Vector) -> float:
    """Get the angle between two vectors.
--- a/propka/version.py
+++ b/propka/version.py
@@ -7,6 +7,7 @@ Contains version-specific methods and parameters.
 TODO - this module unnecessarily confuses the code.  Can we eliminate it?
 """
 import logging
+from propka.atom import Atom
 from propka.hydrogens import setup_bonding_and_protonation, setup_bonding
 from propka.hydrogens import setup_bonding_and_protonation_30_style
 from propka.energy import radial_volume_desolvation, calculate_pair_weight
@@ -98,6 +99,10 @@ class Version:
        """Setup bonding using assigned model."""
        return self.prepare_bonds(self.parameters, molecular_container)

+    def get_hydrogen_bond_parameters(self, atom1: Atom, atom2: Atom) -> tuple:
+        """Get hydrogen bond parameters for two atoms."""
+        raise NotImplementedError("abstract method")
+

 class VersionA(Version):
    """TODO - figure out what this is."""
--- a/tests/test_basic_regression.py
+++ b/tests/test_basic_regression.py
@@ -68,11 +68,11 @@ def run_propka(options, pdb_path, tmp_path):
    """
    options += [str(pdb_path)]
    args = loadOptions(options)
+    cwd = Path.cwd()
    try:
        _LOGGER.warning(
            "Working in tmpdir {0:s} because of PROPKA file output; "
            "need to fix this.".format(str(tmp_path)))
-        cwd = Path.cwd()
        os.chdir(tmp_path)
        parameters = read_parameter_file(args.parameters, Parameters())
        molecule = MolecularContainer(parameters, args)
@@ -148,6 +148,7 @@ def compare_output(pdb, tmp_path, ref_path):
 def test_regression(pdb, options, tmp_path):
    """Basic regression test of PROPKA functionality."""
    path_dict = get_test_dirs()
+    ref_path = None

    for ext in ["json", "dat"]:
        ref_path = path_dict["results"] / f"{pdb}.{ext}"
--- a/tests/test_vector_algebra.py
+++ b/tests/test_vector_algebra.py
@@ -0,0 +1,159 @@
+import propka.vector_algebra as m
+
+import math
+import pytest
+from pytest import approx
+
+RADIANS_90 = math.pi / 2
+
+
+def assert_vector_equal(v1: m.Vector, v2: m.Vector):
+    assert isinstance(v1, m.Vector)
+    assert isinstance(v2, m.Vector)
+    assert v1.x == approx(v2.x)
+    assert v1.y == approx(v2.y)
+    assert v1.z == approx(v2.z)
+
+
+def _matrix4x4_tolist(self: m.Matrix4x4) -> list:
+    return [
+        self.a11, self.a12, self.a13, self.a14, self.a21, self.a22, self.a23, self.a24,
+        self.a31, self.a32, self.a33, self.a34, self.a41, self.a42, self.a43, self.a44
+    ]
+
+def assert_matrix4x4_equal(m1: m.Matrix4x4, m2: m.Matrix4x4):
+    assert isinstance(m1, m.Matrix4x4)
+    assert isinstance(m2, m.Matrix4x4)
+    assert _matrix4x4_tolist(m1) == approx(_matrix4x4_tolist(m2))
+
+
+def test_Vector__init():
+    v = m.Vector()
+    assert v.x == 0.0
+    assert v.y == 0.0
+    assert v.z == 0.0
+    v = m.Vector(12, 34, 56)
+    assert v.x == 12
+    assert v.y == 34
+    assert v.z == 56
+    v1 = m.Vector(atom1=v)
+    assert v1.x == 12
+    assert v1.y == 34
+    assert v1.z == 56
+    v2 = m.Vector(5, 4, 3)
+    v3 = m.Vector(atom1=v2, atom2=v1)
+    assert v3.x == 7
+    assert v3.y == 30
+    assert v3.z == 53
+
+
+def test_Vector__plusminus():
+    v1 = m.Vector(1, 2, 3)
+    v2 = m.Vector(4, 5, 6)
+    v3 = v1 + v2
+    assert v3.x == 5
+    assert v3.y == 7
+    assert v3.z == 9
+    v3 = v1 - v2
+    assert v3.x == -3
+    assert v3.y == -3
+    assert v3.z == -3
+    v4 = -v1
+    assert v4.x == -1
+    assert v4.y == -2
+    assert v4.z == -3
+
+
+def test_Vector__mul__number():
+    v1 = m.Vector(1, 2, 3)
+    assert_vector_equal(v1 * 2, m.Vector(2, 4, 6))
+
+
+def test_Vector__mul__Vector():
+    v1 = m.Vector(1, 2, 3)
+    v2 = m.Vector(4, 5, 6)
+    assert v1 * v2 == 32
+    assert v1.dot(v2) == 32
+    with pytest.raises(TypeError):
+        v1 @ v2  # type: ignore
+
+
+def test_Vector__mul__Matrix4x4():
+    v1 = m.Vector(1, 2, 3)
+    assert_vector_equal(v1 * m.Matrix4x4(), m.Vector())
+    m2 = m.Matrix4x4(0, 1, 0, 0, 20, 0, 0, 0, 0, 0, 300, 0, 0, 0, 0, 1)
+    assert_vector_equal(v1 * m2, m.Vector(2, 20, 900))
+    assert_vector_equal(m2 * v1, m.Vector(2, 20, 900))
+    assert_vector_equal(m2 @ v1, m.Vector(2, 20, 900))
+    with pytest.raises(TypeError):
+        v1 @ m2  # type: ignore
+
+
+def test_Vector__cross():
+    v1 = m.Vector(1, 2, 3)
+    v2 = m.Vector(4, 5, 6)
+    assert_vector_equal(v1**v2, m.Vector(-3, 6, -3))  # TODO deprecate
+    assert_vector_equal(v1.cross(v2), m.Vector(-3, 6, -3))
+    assert_vector_equal(v2.cross(v1), m.Vector(3, -6, 3))
+
+
+def test_Vector__length():
+    v1 = m.Vector(1, 2, 3)
+    assert v1.length() == 14**0.5
+    assert v1.sq_length() == 14
+
+
+def test_Vector__orthogonal():
+    v1 = m.Vector(1, 2, 3)
+    assert v1.dot(v1.orthogonal()) == 0
+
+
+def test_Vector__rescale():
+    v1 = m.Vector(1, 2, 3)
+    v2 = v1.rescale(4)
+    assert v2.length() == 4
+    assert v2.x / v1.x == approx(4 / 14**0.5)
+    assert v2.y / v1.y == approx(4 / 14**0.5)
+    assert v2.z / v1.z == approx(4 / 14**0.5)
+
+
+def test_angle():
+    v1 = m.Vector(0, 0, 1)
+    v2 = m.Vector(0, 2, 0)
+    assert m.angle(v1, v2) == RADIANS_90
+
+
+def test_angle_degrees():
+    v1 = m.Vector(0, 0, 3)
+    v2 = m.Vector(5, 0, 0)
+    assert m.angle_degrees(v1, v2) == 90
+
+
+def test_signed_angle_around_axis():
+    v1 = m.Vector(0, 0, 3)
+    v2 = m.Vector(5, 0, 0)
+    v3 = m.Vector(0, 1, 0)
+    assert m.signed_angle_around_axis(v1, v2, v3) == -RADIANS_90
+    v1 = m.Vector(0, 2, 3)
+    v2 = m.Vector(5, 4, 0)
+    assert m.signed_angle_around_axis(v1, v2, v3) == -RADIANS_90
+
+
+def test_rotate_vector_around_an_axis():
+    v1 = m.Vector(0, 0, 3)
+    v2 = m.Vector(3, 0, 0)
+    v3 = m.Vector(0, -1, 0)
+    v4 = m.rotate_vector_around_an_axis(RADIANS_90, v3, v2)
+    assert_vector_equal(v4, v1)
+
+
+def test_rotate_atoms_around_z_axis():
+    m_rot = m.rotate_atoms_around_z_axis(-RADIANS_90)
+    assert_matrix4x4_equal(m_rot,
+                           m.Matrix4x4(0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1))
+
+
+def test_rotate_atoms_around_y_axis():
+    m_rot = m.rotate_atoms_around_y_axis(RADIANS_90)
+    assert_matrix4x4_equal(m_rot,
+                           m.Matrix4x4(0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1))