dynamics 0.2.0

Molecular dynamics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

//! Contains setup code, including applying forcefield data to our specific
//! atoms.

// Notes to square away the 3 "atom name" / "Amber atom type" / "force field type" keys.
// This guide shows Type 1. https://emleddin.github.io/comp-chem-website/AMBERguide-AMBER-atom-types.html,
//
// Update: "Type 1" = "type_in_res" in our code now. "Type 2" = "ff_type" for AAs, and "Type 3" = "ff_type" for small mols.
//
// Type 1 Examples: "CA", "HA", "CZ", "HB3", "HH22", HZ2", "N", "H", "HG3", "O", "CD", "C", "HG23", "CG", "CB", "CG1", "HE2", "HB3",
// Type 1 Sources: `amino19.lib`, col 0. mmCIF atom coordinate files.
//
// Type 2 Examples:  "HC", "C8", "HC", "H"(both), "XC", "N"(both), "H"(both), "H1", "CT", "OH", "HO", "2C",
// Type 2 Sources: `amino19.lib` (AA/protein partial charges), col 1. `frcmod.ff19SB`. (AA/protein params)
//
// Small Mol/lig:
// Type 3 Examples: "oh", "h1", "ca", "o", "os", "c6", "n3", "c3"
// Type 3 Sources: `.mol2` generated by Amber. (Small mol coordinates and partial charges) `gaff2.dat` (Small molg params)
//
// MOl2 for ligands also have "C10", "O7" etc, which is ambiguous here, and not required, as their params
// use Type 3, which is present. It's clear what to do for ligand
//
// Best guess: Type 1 identifies labels within the residue only. Type 2 (AA) and Type 3 (small mol) are the FF types.

use std::{collections::HashSet, fmt};

#[cfg(feature = "encode")]
use bincode::{Decode, Encode};
use bio_files::{
    gromacs::mdp::{ConstraintAlgorithm, Constraints},
    md_params::{
        AngleBendingParams, BondStretchingParams, ForceFieldParams, ForceFieldParamsIndexed,
        LjParams, MassParams,
    },
};
use itertools::Itertools;
use na_seq::Element::Hydrogen;

use crate::{AtomDynamics, MdState, ParamError};

/// Add items from one parameter set to the other. If there are duplicates, the second set's overrides
/// the baseline.
pub fn merge_params(baseline: &ForceFieldParams, add_this: &ForceFieldParams) -> ForceFieldParams {
    let mut merged = baseline.clone();

    merged.mass.extend(add_this.mass.clone());
    merged.lennard_jones.extend(add_this.lennard_jones.clone());

    merged.bond.extend(add_this.bond.clone());
    merged.angle.extend(add_this.angle.clone());
    merged.dihedral.extend(add_this.dihedral.clone());
    merged.improper.extend(add_this.improper.clone());

    merged
}

/// We use this variant in the configuration API. Deferrs to `HydrogenConstraintInner` for holding
/// constraints.
#[cfg_attr(feature = "encode", derive(Encode, Decode))]
#[derive(Clone, Copy, PartialEq, Debug)]
pub enum HydrogenConstraint {
    /// A linear constraint solver. Similar to GROMACS' LINCS. This is slightly faster,
    /// and slightly more stable than SHAKE.
    Linear {
        order: u8,
        iter: u8,
    },
    Shake {
        shake_tolerance: f32,
    },
    /// Uses Shake and Rattle to fix the hydrogen positions. This allows for a larger timestep,
    /// e.g. 2fs instead of 1fs.
    /// Uses the same bonded parameters as elsewhere: A spring model
    Flexible,
}

impl Default for HydrogenConstraint {
    fn default() -> Self {
        // HydrogenConstraint::Shake {
        //     shake_tolerance: SHAKE_TOL_DEFAULT,
        // }
        // GROMACS defaults here.
        HydrogenConstraint::Linear { order: 4, iter: 1 }
    }
}

impl HydrogenConstraint {
    pub fn to_gromacs(&self) -> Constraints {
        match self {
            Self::Linear { order, iter } => Constraints::HBonds(ConstraintAlgorithm::Lincs {
                order: *order,
                iter: *iter,
            }),
            Self::Shake {
                shake_tolerance: tolerance,
            } => Constraints::HBonds(ConstraintAlgorithm::Shake { tol: *tolerance }),
            Self::Flexible => Constraints::None,
        }
    }
}

impl fmt::Display for HydrogenConstraint {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            HydrogenConstraint::Linear { .. } => {
                write!(f, "Constrained (LINCS)")
            }
            HydrogenConstraint::Shake { shake_tolerance: _ } => {
                write!(f, "Constrained (SHAKE)")
            }
            HydrogenConstraint::Flexible => write!(f, "Flexible"),
        }
    }
}

impl MdState {
    /// We use this to set up optimizations defined in the Amber reference manual. `excluded` deals
    /// with sections were we skip coulomb and Vdw interactions for atoms separated by 1 or 2 bonds. `scaled14` applies a force
    /// scaler for these interactions, when separated by 3 bonds.
    pub(crate) fn setup_nonbonded_exclusion_scale_flags(&mut self) {
        // Helper to store pairs in canonical (low,high) order
        let push = |set: &mut HashSet<(usize, usize)>, i: usize, j: usize| {
            if i < j {
                set.insert((i, j));
            } else {
                set.insert((j, i));
            }
        };

        // 1-2
        for indices in &self.force_field_params.bonds_topology {
            push(&mut self.pairs_excluded_12_13, indices.0, indices.1);
        }

        // 1-3
        for indices in self.force_field_params.angle.keys() {
            push(&mut self.pairs_excluded_12_13, indices.0, indices.2);
        }

        // 1-4. We do not count improper dihedrals here.
        for indices in self.force_field_params.dihedral.keys() {
            push(&mut self.pairs_14_scaled, indices.0, indices.3);
        }

        // Make sure no 1-4 pair is also in the excluded set
        for p in &self.pairs_14_scaled {
            self.pairs_excluded_12_13.remove(p);
        }
    }
}