Description:

A critical bug has been identified in the energy calculation module where the hydrogen bond energy term (E_hbond) is consistently calculated as zero or is incorrect. This significantly impacts the accuracy of the total system energy, leading to suboptimal side-chain placements.

The root cause of this issue stems from a multi-stage failure in the hydrogen bond parameterization and scoring pipeline:

1. Forcefield Parsing: The current parser for forcefield files (.toml) correctly reads the [hbond] table, but it fails to interpret the Donor-Acceptor relationships from the keys (e.g., "O_3-O_3M"). As a result, the application has no internal knowledge of which atom types can act as hydrogen bond donors or acceptors.

2. Atom Parameterization: The Parameterizer is responsible for assigning physicochemical properties to each atom. Due to the lack of donor/acceptor information from the parser, it cannot correctly assign a hydrogen bond role (hbond_type_id) to atoms. The current logic is flawed, resulting in all atoms being assigned a default value of -1 (not involved in H-bonding).

3. Energy Scoring: The Scorer relies on the hbond_type_id to identify potential hydrogen bond interactions (D-H...A). Since this ID is never set correctly, the scoring function fails to find any valid hydrogen bond triplets and returns an energy of zero.

Tasks:

• Phase 1: Enhance Forcefield Parameter Loading (scream-core/src/core/forcefield/params.rs)

  • In the NonBondedParams struct, add two new HashSet<String> fields: hbond_donors and hbond_acceptors. These should be marked with #[serde(skip)] as they will be populated manually.
  • Modify the Forcefield::load_non_bonded function. After successfully parsing the TOML file, iterate through the keys of the hbond HashMap. For each key (e.g., "D-A"), parse the donor and acceptor types and populate the hbond_donors and hbond_acceptors sets accordingly.

• Phase 2: Implement Correct H-Bond Role Parameterization (scream-core/src/core/forcefield/parameterization.rs)

  • Create a new private method assign_hbond_roles(&self, system: &mut MolecularSystem) within the Parameterizer.
  • In assign_hbond_roles, implement the logic to iterate through all atoms in the system and assign their hbond_type_id based on the following rules:
    • Acceptor (ID=1): If an atom's force_field_type is present in the hbond_acceptors set.
    • Donor Hydrogen (ID=0): If an atom's name starts with 'H', it has exactly one covalent bond, and its bonded neighbor's force_field_type is present in the hbond_donors set.
    • None (ID=-1): For all other atoms.
  • Call the new assign_hbond_roles method at the end of parameterize_system to ensure it runs after initial atom roles and types are set.
  • Update parameterize_rotamer to apply the same H-bond role assignment logic for individual rotamers.
  • Remove the old, incorrect hbond_type_id assignment logic.

• Phase 3: Refactor H-Bond Energy Calculation (scream-core/src/core/forcefield/scoring.rs)

  • Rewrite the Scorer::calculate_hbond_one_way function to use the newly assigned hbond_type_id.
  • The logic should first identify all potential donor hydrogens (hbond_type_id == 0) in one group.
  • For each donor hydrogen (H), use system.get_bonded_neighbors to find its covalently attached donor heavy atom (D).
  • Then, iterate through the second group to find all potential acceptor atoms (hbond_type_id == 1).
  • For each valid D-H...A triplet found, perform the necessary geometric checks (distance, angle).
  • Look up the correct HBondParam from the forcefield using the force_field_type of the donor (D) and acceptor (A).
  • Call EnergyCalculator::calculate_hbond with the correct atoms and parameters to compute the energy contribution.