Entire namespace¶

namespace mdk¶

Typedefs

using Scalars = Eigen::VectorXd¶: A list of scalars (i.e. doubles).

typedef Eigen::Vector3d Vector¶: A vector (writing Eigen::Vector3d everywhere is annoying).

typedef Vector const &VRef¶: Constant reference to a vector. Given the size of a Vector, passing copies wouldn’t be inefficient; it’s here chiefly to hide the warnings.

using Bytes = std::vector<int8_t>¶: List of small integers (usually we use this to store various enums and true/false values).

using Integers = std::vector<int>¶

using Pairs = std::vector<std::pair<int, int>>¶: List of pairs (used in storing VL pairs).

using VectorBase = Eigen::Matrix<double, 3, Eigen::Dynamic, Eigen::ColMajor>¶

Enums

enum AminoAcidIdx¶

An underlying amino acid index.

Values:

enumerator ALA¶

enumerator ARG¶

enumerator ASN¶

enumerator ASP¶

enumerator CYS¶

enumerator GLU¶

enumerator GLN¶

enumerator GLY¶

enumerator HIS¶

enumerator ILE¶

enumerator LEU¶

enumerator LYS¶

enumerator MET¶

enumerator PHE¶

enumerator PRO¶

enumerator SER¶

enumerator THR¶

enumerator TRP¶

enumerator TYR¶

enumerator VAL¶

enum ContactTypeIdx¶

An underlying contact type index.

Values:

enumerator NAT¶

enumerator NAT_BB¶

enumerator NAT_BS¶

enumerator NAT_SB¶

enumerator NAT_SS¶

enumerator SSBOND¶

enum PairType¶

A “pair type”, as distinguished by a number of force fields, such as HeuresticBondAngles or HeuresticDihedralAngles.

Values:

enumerator GG¶

enumerator GP¶

enumerator GX¶

enumerator PG¶

enumerator PP¶

enumerator PX¶

enumerator XG¶

enumerator XP¶

enumerator XX¶

enum ResTypeIdx¶

Underlying residue type index.

Values:

enumerator ALA¶

enumerator ARG¶

enumerator ASN¶

enumerator ASP¶

enumerator CYS¶

enumerator GLU¶

enumerator GLN¶

enumerator GLY¶

enumerator HIS¶

enumerator ILE¶

enumerator LEU¶

enumerator LYS¶

enumerator MET¶

enumerator PHE¶

enumerator PRO¶

enumerator SER¶

enumerator THR¶

enumerator TRP¶

enumerator TYR¶

enumerator VAL¶

Functions

std::vector<PairType> pairTypes()¶

PairType pairType(AminoAcid const &acid1, AminoAcid const &acid2)¶

inline double angle(VRef v1, VRef v2, VRef v3)¶

inline double dihedral(VRef v1, VRef v2, VRef v3, VRef v4)¶

double asphericity(Vectors const &v)¶

std::string_view view(std::string const &s, int i, int j)¶

char view(std::string const &s, int i)¶

std::string_view ltrim(std::string_view s)¶

std::string_view rtrim(std::string_view s)¶

std::string_view trim(std::string_view s)¶

std::string line(std::istream &is)¶

void skipLine(std::istream &is)¶

std::stringstream lineStream(std::istream &is)¶

Variables

constexpr int numOfPTs = 9¶

Unit f77unit = 1.0¶

Unit angstrom = f77unit / 5.0¶

Unit nanometer = 10.0 * angstrom ¶

Unit meter = nanometer * 1.0e9¶

Unit nanosecond = 1.0¶

Unit tau = 1.0 * nanosecond ¶

Unit microsecond = nanosecond * 1.0e3¶

Unit millisecond = nanosecond * 1.0e6¶

Unit second = nanosecond * 1.0e9¶

Unit atom = 1.0¶

Unit mol = 6.02214076e23 / atom ¶

Unit eps = 1.0¶

Unit kcal = eps * mol / 1.57824959¶

Unit Joule = kcal / 4184.0¶

Unit eps_kB = 1.0¶

Unit kB = eps / eps_kB ¶

Unit Kelvin = 1.380649e-23 * Joule / kB ¶

Unit kg = Joule * second * second / (meter * meter)¶

Unit au = kg * 0.99999999965e-3 / mol ¶

Unit f77mass = eps * tau * tau / (f77unit * f77unit)¶: In the Fortran version of the model, distance of f77unit, time of tau, energy of eps and the average mass of an aminoacid were units; these are however incongruent, it’s a confirmed bug.

Unit echarge = 1.0¶

Unit Coulomb = echarge / 1.602176634e-19¶

Unit Ampere = Coulomb / second ¶

Unit cspeed = 299792458.0 * meter / second ¶

Unit Henry = kg * meter * meter / (second * second * Ampere * Ampere)¶

Unit mu_0 = 1.25663706212e-6 * Henry / meter ¶

Unit epsilon_0 = 1.0 / (mu_0 * cspeed * cspeed)¶

Unit rad = 1.0¶

Unit degree = (2.0 * M_PI / 360.0) * rad ¶

class Chains¶

#include <Chains.hpp>

This data structure contains all the data pertaining to chains, especially the pairs, triples etc. in the model. It’s here due to similarity of computation and some degree of reuse of the same data.

The definitions of pairs, triples etc. follow convention from CPC14.pdf.

Public Functions

Chains() = default¶

explicit Chains(Model const &model)¶

Generate a Chains structure for a given model.

Parameters: model – The model, chain data whereof we wish to obtain.

inline bool sepByAtLeastN(int i1, int i2, int n) const¶

Check bond distance between residues.

Parameters

i1 – First residue
i2 – Second residue
n – Bond offset

Returns

true if residues i1 and i2 are either in separate chains or in a single chain but separated by at least n bonds.

Bytes tuples(int k) const¶

Create an array of tuple indicators.

Parameters: k – Tuple size
Returns: An array of size model->n, where tuples(k)[i] is 1 if a tuple (i-k+2, …, i+1) is entirely contained in a chain, and 0 otherwise.

Bytes nativeTuples(int k) const¶

Create an array of native tuple indicators.

Parameters: k – Tuple size
Returns: An array of size model->n where tuples(k)[i] is 1 if a tuple (i-k+2, …, i+1) is entirely contained in a structured part, and 0 otherwise.

Public Members

Model const *model = nullptr¶: Const ptr to an underlying model (necessarily it must have at least as long a lifetime as this class).

Integers chainIdx¶: chainIdx[i] is the index of the chain for residue i, or -1 if a residue is in no chain.

Bytes isTerminal¶: isTerminal[i] = 1 if the i’th residue is first or last in the chain (and thus lacks a neighbor), or 0 otherwise.

Bytes pairs¶: pairs[i] = 1 if the pair (i, i+1) is entirely in one chain, or 0 otherwise. Result of tuples(2). Used in harmonic tether forcefield.

Bytes triples¶: triples[i] = 1 if the triple (i-1, i, i+1) is entirely in one chain, or 0 otherwise. Result of tuples(3). Used for QA (in computing n, h) and for bond angle potentials.

Bytes quads¶: quads[i] = 1 if the quadruple (i-2, i-1, i, i+1) is entirely in one chain, or 0 otherwise. Result of tuples(4). Used in dihedral potential.

Bytes nativeTriples¶: nativeTriples[i] = 1 if the triple (i-1, i, i+1) is entirely in a structured part, or 0 otherwise. Result of nativeTuples(3). Used in native bond angle potential.

Bytes nativeQuads¶: nativeQuads[i] = 1 if the quadruple (i-2, i-1, i, i+1) is entirely in a structured part, or 0 otherwise. Result of nativeTuples(4). Used in native dihedral potentials.

class Charges : public Scalars ¶

#include <Charges.hpp>

A data structure containing the charges of residues in a model, as derived from a parameter file, in an array form.

Public Functions

Charges() = default¶

Charges(Model const &model, param::Parameters const &params)¶

class DataFactory¶

#include <DataFactory.hpp>

A generic provider of data classes derived from a model and a parameter file. This class is provided chiefly for convenience.

Public Functions

inline DataFactory(Model const &model, param::Parameters const &params)¶

template<typename Data> inline Data const &data()¶

Extract a data class of type Data, if such is defined; otherwise, one will most likely get a linker error, as relevant template specialzations of /p create are placed in a .cpp file.

Template Parameters: Data – Type of data class to retrieve
Returns: A const reference to a stored data class.

Private Functions

template<typename Data> Data create() const¶

Underlying data class constructor.

Template Parameters: Data – Data class to construct
Returns: Constructed data class.

template<> Chains create() const¶

template<> Charges create() const¶

template<> Masses create() const¶

template<> Types create() const¶

template<> Model create() const¶

Private Members

Model const *model¶

param::Parameters const *params¶

std::unordered_map<std::type_index, std::shared_ptr<void>> savedData¶: A realization of a heterogeneous container. It holds values based on the types of elements, i.e. at most one value of each type can be stored. Using std::shared_ptr<void> allows for safe destruction of stored entities.

class Masses : public Scalars ¶

#include <Masses.hpp>

Data class containing masses associated with residues, in an array form.

Public Functions

Masses() = default¶

explicit Masses(Model const &model)¶

class Types : public Eigen::Matrix<ResType, Eigen::Dynamic, 1>¶

#include <Types.hpp>

Data class containing types of residues, in an array form.

Public Functions

Types() = default¶

explicit Types(Model const &model)¶

class Vectors : public VectorBase ¶

#include <Vectors.hpp>

This class is a sort of a decorator over standard Eigen::Matrix3Xd. The chief difference is that this class offers an interface more appropriate for dealing with it as though it were a list of vectors (for example, one needs not use m.col(i) to get t’th vector in the matrix. Also, now size() corresponds to the number of vectors and not the elements of the array. The array itself is not (easily) resizeable.

Public Functions

Vectors() = default¶

inline explicit Vectors(int n)¶

Create an unitialized list of n vectors.

Parameters: n – Number of vectors

inline Vectors(int n, Vector const &init)¶

Create an initialized list of n vectors.

Parameters

n – Number of vectors
init – Initial value

inline auto vectorwise()¶

Returns an iterator over constituent vectors.

Returns: An interator over constituent vectors.

inline auto vectorwise() const¶

Returns a const iterator over constituent vectors.

Returns: A const interator over constituent vectors.

inline int size() const¶

Returns: A number of vectors.

inline auto operator[](int i)¶

Access to vector.

Parameters: i – Index of a vector to access
Returns: A slice of a matrix corresponding to i’th vector.

inline auto operator[](int i) const¶

Const access to vector.

Parameters: i – Index of a vector to access
Returns: A const slice of a matrix corresponding to i’th vector.

class BondAngles : public mdk::Force ¶

#include <BondAngles.hpp>

A force field for bond angles (both heurestic and native variants, depending on whether a triple has an associated native bond angle (in such cases the native variant supercedes the heurestic variant).

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force field to a simulation.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

“Asynchronous” part of the potential. Here it’s the only part.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Private Members

HeuresticBA *heurBA = nullptr¶: The “heurestic” part of the potential.

NativeBA *natBA = nullptr¶: The “native” part of the potential.

Bytes inRange¶: Whether a triple (i-1, i, i+1) is connected, i.e. in one chain.

Friends

friend class HeuresticBA

friend class NativeBA

class HeuresticBA : public mdk::SimulVar ¶

#include <HeuresticBA.hpp>

“Heurestic” part of the bond angle potential, applied when a triple has no associated native bond angle.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the part to the simulation - in particular, find (and possibly create) BondAngles force field and add itself to it.

Parameters: simulation – Simulation to bind to.

inline void term(int i, double theta, double &V, double &dV_dth) const¶

Term of the formula for bond angle potential. Note: it’s inline in order for the compiler to inline it in BondAngles.cpp file.

Parameters

i – The index i in the triple (i-1, i, i+1).
theta – Value of the bond angle between i-1, i and i+1.
V – Potential energy reference to add to.
dV_dth – Derivative of potential energy wrt the angle theta to add to.

Private Members

double coeff[numOfPTs][D + 1]¶: Polynomial coefficients.

Bytes angleTypes¶: List of “types of triples”, or rather pair types for (i, i+1) for a triple (i-1, i, i+1).

Private Static Attributes

static constexpr const int D = 6¶: Degree of the polynomial.

Friends

friend class BondAngles

class NativeBA : public mdk::SimulVar ¶

#include <NativeBA.hpp>

“Native” part of the bond angle potential, applied when a triple has an associated native bond angle.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the part to the simulation - in particular, find (and possibly create) BondAngles force field and add itself to it.

Parameters: simulation – Simulation to bind to.

inline void term(int i, double theta, double &V, double &dV_dth) const¶

Term of the formula for bond angle potential. Note: it’s inline in order for the compiler to inline it in BondAngles.cpp file.

Parameters

i – The index i in the triple (i-1, i, i+1).
theta – Value of the bond angle between i-1, i and i+1.
V – Potential energy reference to add to.
dV_dth – Derivative of potential energy wrt the angle theta to add to.

Private Members

Bytes isNative¶: isNative[i] is 1 if (i-1, i, i+1) has an associated native bond angle, 0 otherwise.

Scalars theta0¶: theta0[i] contains the native bond angle whenever such is defined. For triples (i-1, i, i+1) which have no such associated angle, the value is undefined.

double CBA = 30.0 * eps / pow(rad, 2)¶: Value of $k_\theta$.

Friends

friend class BondAngles

class Chirality : public mdk::Force ¶

#include <Chirality.hpp>

Chirality potential for the models with a native structure (in particular the positions of the residues). It’s supposed to be a more replacement for bond and dihedral potentials, if the native structure is provided.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force field to the simulation.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Public Members

double e_chi = 1.0 * eps ¶: The amplitude of the potential.

Private Members

Scalars d0_cube_inv¶: A list of cube inverses of $d_0$, where $d_0$ is the length of $v_i = r_{i+1} - r_i$ in the native structure.

Scalars C_nat¶: A list of the values of C_i in the native structure.

Bytes inRange¶: inRange[i] = 1 if a quadruple (i-2, i+1, i, i+1) is connected, i.e. in a single chain.

class ComplexNativeDihedral : public mdk::NativeDihedralBase ¶

#include <ComplexNativeDihedral.hpp>

The complex variant of the native part of the dihedral angle potential.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the part to the simulation - in particular, find (and possibly create) DihedralAngles force field and add itself to it.

Parameters: simulation – Simulation to bind to.

inline void term(int i, double phi, double &V, double &dV_dphi) const¶

Term of the formula for the dihedral angle potential of a single quadruple. Note: it’s inline in order for the compiler to inline it in DihedralAngles.cpp file.

Parameters

i – The index i in the quadruple (i-2, i-1, i, i+1).
phi – Dihedral angle between planes defined by sequences i-2, i-1, i and i-1, i, i+1.
V – Potential energy reference to add to.
dV_dphi – Derivative of potential energy wrt the angle phi to add to.

Public Members

double CDA = 0.66 * eps / pow(rad, 2.0)¶

double CDB = 0.66 * eps / pow(rad, 2.0)¶

Friends

friend class DihedralAngles

class DihedralAngles : public mdk::Force ¶

#include <DihedralAngles.hpp>

A force field for dihedral angles, composed of:

the “native” part, which exists in two mutually exclusive variants (simple and complex);
the “heurestic” part. The native part supercedes heurestic part whenever native dihedral angle is defined for a quadruple.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force field to a simulation.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

“Asynchronous” part of the potential. Here it’s the only part.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Private Members

std::variant<std::monostate, ComplexNativeDihedral*, SimpleNativeDihedral*> natDih = std::monostate()¶: The “native” part, i.e. either nothing, complex variant or the simple variant.

HeuresticDihedral *heurDih = nullptr¶: The “heurestic” part of the potential.

Bytes inRange¶

Friends

friend class ComplexNativeDihedral

friend class SimpleNativeDihedral

friend class HeuresticDihedral

class HeuresticDihedral : public mdk::SimulVar ¶

#include <HeuresticDihedral.hpp>

Heurestic part of the dihedral potential, applied when with the quadruple (i-2, i-1, i, i+1) is not associated a native dihedral angle.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the part to the simulation - in particular, find (and possibly create) DihedralAngles force field and add itself to it.

Parameters: simulation – Simulation to bind to.

inline void term(int i, double phi, double &V, double &dV_dphi) const¶

Term of the formula for the dihedral angle potential of a single quadruple. Note: it’s inline in order for the compiler to inline it in DihedralAngles.cpp file.

Parameters

i – The index i in the quadruple (i-2, i-1, i, i+1).
phi – Dihedral angle between planes defined by sequences i-2, i-1, i and i-1, i, i+1.
V – Potential energy reference to add to.
dV_dphi – Derivative of potential energy wrt the angle phi to add to.

Private Members

double coeff[numOfPTs][6]¶: Coefficients per pair type (where pair is (i-1, i) from the quadruple.

Eigen::Matrix<int8_t, Eigen::Dynamic, 1> angleTypes¶: Types of angles (i.e. pair types for each pair (i-1, i)).

Friends

friend class DihedralAngles

class NativeDihedralBase : public mdk::SimulVar ¶

#include <NativeDihedralBase.hpp>

A common part of simple and complex native dihedral variants. In particular, it stores whether a quadruple (i-2, i-1, i, i+1) has an associated native dihedral angle, and its value if it is the case.

Subclassed by mdk::ComplexNativeDihedral, mdk::SimpleNativeDihedral

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the variant to the simulation. Technically one shouldn’t add this base to the simulation itself; it serves as a prototype for bind in ComplexNativeDihedral and SimpleNativeDihedral in particular, it retrieves and computes isNative and phi0.

Parameters: simulation –

Protected Attributes

Bytes isNative¶: isNative[i] = 1 if a quadruple (i-2, i-1, i, i+1) has an associated native dihedral angle, 0 otherwise.

Scalars phi0¶: phi0[i] has the value of an associated native dihedral angle, if such exists for the quadruple (i-2, i-1, i, i+1); otherwise, the value is undefined.

class SimpleNativeDihedral : public mdk::NativeDihedralBase ¶

#include <SimpleNativeDihedral.hpp>

The simple variant of the native part of the dihedral angle potential.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the part to the simulation - in particular, find (and possibly create) DihedralAngles force field and add itself to it.

Parameters: simulation – Simulation to bind to.

inline void term(int i, double phi, double &V, double &dV_dphi) const¶

Term of the formula for the dihedral angle potential of a single quadruple. Note: it’s inline in order for the compiler to inline it in DihedralAngles.cpp file.

Parameters

i – The index i in the quadruple (i-2, i-1, i, i+1).
phi – Dihedral angle between planes defined by sequences i-2, i-1, i and i-1, i, i+1.
V – Potential energy reference to add to.
dV_dphi – Derivative of potential energy wrt the angle phi to add to.

Public Members

double CDH = 3.33 * eps / pow(rad, 2)¶

Friends

friend class DihedralAngles

class ConstDH : public mdk::ESBase ¶

#include <ConstDH.hpp>

Debye-Hueckel screeened electrostatic potential with constant electric permittivity of the medium.

Public Functions

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Public Members

double screeningDist = 10.0 * angstrom ¶

double permittivity = 80.0 * epsilon_0 ¶

Protected Functions

virtual vl::Spec spec() const override¶

Generate a VL spec.

Returns: Generated VL spec.

class ESBase : public mdk::NonlocalForce ¶

#include <ESBase.hpp>

A base class for the Debye-Hueckel screened electrostatic potentials.

Subclassed by mdk::ConstDH, mdk::RelativeDH

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the class to the simulation. It initializes charge and adds the force to the VL list. Note: this class (ESBase) shouldn’t be added to the simulation class, only the derived classes.

Parameters: simulation – Simulation to bind to.

virtual void vlUpdateHook() override¶: An action to be performed when the Verlet list is updated. Here we want to recreate pairs, i.e. filter the Verlet list so that only charged pairs are present.

Protected Attributes

std::vector<Contact> pairs¶: A local, filtered Verlet list of pairs. Since most residues are not charged, and pairs even moreso, it benefits us to keep a local list and not check the pairs which we would know cannot interact with each other.

Eigen::Matrix<int8_t, Eigen::Dynamic, 1> charge¶: A list of charges of the residues (in units of echarge).

struct Contact¶

#include <ESBase.hpp>

A struct detailing a contact between pairs. We store the indices and the product of the charges.

Public Members

int i1¶: Index of the first residue.

int i2¶: Index of the second residue.

double q1_x_q2¶: Product of residues’ charges.

class RelativeDH : public mdk::ESBase ¶

#include <RelativeDH.hpp>

Debye-Hueckel screened electrostatic potential with a distance-dependent electric permittivity (specifically eps_r = r_0/r).

Public Functions

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Public Members

double screeningDist = 10.0 * angstrom ¶

double r0 = 4.0 * angstrom ¶

Protected Functions

virtual vl::Spec spec() const override¶

Generate a VL spec.

Returns: Generated VL spec.

class Force : public mdk::SimulVar ¶

#include <Force.hpp>

An abstract interface for the force fields. Certain forces (QA in particular) have a multi-stage evaluation. Thus we separated the computation to an “asynchronous” and “synchronous” parts, where the latter get executed in order of addition to the Simulation object.

Subclassed by mdk::BondAngles, mdk::Chirality, mdk::DihedralAngles, mdk::NonlocalForce, mdk::SolidWall, mdk::Tether

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force to the simulation. This class shouldn’t be actually added to the simulation, but rather serves as a prototype for actual forces; in particular it saved state from the Simulation object.

Parameters: simulation –

virtual void asyncPart(Dynamics &dynamics)¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

virtual void syncPart(Dynamics &dynamics)¶

Synchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Protected Attributes

State const *state = nullptr¶: A const pointer to the State of simulation. The const-ness is a declaration that the State is only modified by the integrator, and thus asyncPart can be safely executed.

class NativeContacts : public mdk::NonlocalForce ¶

#include <NativeContacts.hpp>

Go model potential. It is somewhat special as it interferes with the Verlet list and with the quasi-adiabatic potential. It should go first in the list of nonlocal forces added to the simulation object.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force to the simulation, and also to the Verlet list.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

virtual void vlUpdateHook() override¶: Action to perform when the Verlet list is updated. We want to (a) retrieve the pairs that are in native contact into a list (curPairs), (b) remove those pairs from the global Verlet list.

Private Functions

virtual vl::Spec spec() const override¶

Generate a spec for the Verlet list.

Returns: Generated spec.

Private Members

Pairs newVL¶: A “copy” of the global Verlet list to be swapped with it when we update it; we want to do the swap so as to avoid excessive allocation of memory when the program runs.

std::vector<Contact> allContacts¶: List of all native contacts.

std::vector<Contact> curPairs¶: List of native contacts that are within the cutoff distance.

double depth = 1.0 * eps ¶: Depth of the Lennard-Jones potential, with which the natively connected residues interact.

struct Contact¶

#include <NativeContacts.hpp>

A struct with contact data.

Public Members

int i1¶: First residue.

int i2¶: Second residue.

double r_min¶: Lennard-Jones potential minimal distance.

class NonlocalForce : public mdk::Force ¶

#include <NonlocalForce.hpp>

A nonlocal force interface. It differs from the standard force interface in that it must provide some specs for the VL list, and needs to implement a hook to run when a VL list is updated.

Subclassed by mdk::ESBase, mdk::NativeContacts, mdk::PauliExclusion, mdk::PseudoImproperDihedral, mdk::QuasiAdiabatic

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the nonlocal force to the simulation object. This class shouldn’t be bound to the simulation object, rather invoked by derived classes. In particular, it fetches Verlet list into vl.

Parameters: simulation – Simulation to bind to.

virtual void vlUpdateHook() = 0¶: An action to be performed when Verlet list is updated; usually one copies the relevant pairs into a local list, perhaps with some modifications. The actions are invoked in the order of being added to the Verlet list.

Protected Functions

virtual vl::Spec spec() const = 0¶

Compute the spec to register to the VL list.

Returns: Generated VL spec.

void installIntoVL()¶: Registers the spec into the Verlet list. Note: it must be executed only one all the data required for its generation has been fetched, i.e. usually at the end of bind.

Protected Attributes

mutable vl::Spec savedSpec¶: A copy of the generated spec(), for use in the computation of forces (for example it holds the cutoff distance).

vl::List *vl = nullptr¶: A pointer to the VL list. It’s not const chiefly because the update hooks may modify the list, in particular create exclusions, like for example the native contacts’ forcefield removing the native contacts

class PauliExclusion : public mdk::NonlocalForce ¶

#include <PauliExclusion.hpp>

A “Pauli exclusion force” that keeps the residues from overlapping.

Public Functions

PauliExclusion()¶

virtual void bind(Simulation &simulation) override¶

Bind the object to the simulation.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

virtual void vlUpdateHook() override¶: An action performed when a Verlet list is reconstructed; here we simply copy the list verbatim.

Public Members

ShiftedTruncatedLJ stlj¶: A shifted and truncated version of the Lennard-Jones potential that is used as the actual force.

Protected Functions

virtual vl::Spec spec() const override¶

Generate a VL spec.

Returns: Generated VL spec.

Private Members

Pairs exclPairs¶: A local copy of the Verlet list.

class LambdaPeak¶

#include <PseudoImproperDihedral.hpp>

An object representing either of the three lambda functions. It is theoretically a Gaussian function, but we implement it here in terms of computationally faster functions (either cosine or an “algebraic version” thereof, with a cutoff at the far ends).

Public Functions

bool supp(double psi) const¶

Check whether a value $psi$ is in the support of the function.

Parameters: psi – Value to test.
Returns: true if $lambda(psi) != 0$, false otherwise.

void eval(double psi, double &L, double &dL_dpsi) const¶

Compute the value of the function.

Parameters

psi – Argument, i.e. the improper dihedral angle
L – Value of the lambda function at psi,
dL_dpsi – Derivative of L wrt psi

Public Members

double alpha¶: The “width” of the lambda function.

double psi0¶: The “center” of the lambda function.

bool cosineVersion¶: Whether we should use the cosine version of the function.

class PseudoImproperDihedral : public mdk::NonlocalForce ¶

#include <PseudoImproperDihedral.hpp>

The “pseudo-improper-dihedral” custom potential. A more computationally expensive but more accurate and Hamiltonian alternative to the quasi-adiabatic potential.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the object to the simulation. Here we also fetch from it the types and chains.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

virtual void vlUpdateHook() override¶: An action to be performed when the Verlet list is updated. We essentially copy the residues but exclude the pairs within bond distance of less than 4, or at the ends of their chains.

Public Members

LambdaPeak bb_pos¶: Lambda function corresponding to the bb+ part of the potential.

LennardJones bb_pos_lj¶: Lennard-Jones potential corresponding to the bb+ part of the potential.

LambdaPeak bb_neg¶: Lambda function corresponding to the bb- part of the potential.

LennardJones bb_neg_lj¶: Lennard-Jones potential corresponding to the bb- part of the potential.

LambdaPeak ss¶: Lambda function corresponding to the ss part of the potential.

SidechainLJ ss_ljs[AminoAcid::N][AminoAcid::N]¶: An array of sidechain L-J potentials for each pair of types of amino acids interacting.

Protected Functions

virtual vl::Spec spec() const override¶

Generate spec for the Verlet list.

Returns: Generated spec.

Private Functions

void deriveAngles(vl::PairInfo const &pair, double psi[2], Vector dpsi_dr[2][6]) const¶

Derive the angle data (specifically the angles and derivatives) pertaining to the pair.

Parameters

pair – Pair of residues.
psi – An array, where psi[0] is the improper dihedral angle between i_1-1, i_1, i_1+1 and i_2, and psi[1] is the improper dihedral angle between i_2-1, i_2, i_2+1 and i_1.
dpsi_dr – An array, where dpsi_dr[0][0..2] are the derivatives d psi_1/d r_{i_1-1 .. i_1+1} and dspi_dr[0][3..6] are the derivatives d psi_1/d r_{i_2-1 .. i_2+1}; dpsi_dr[1][0..3] is defined as for 0, but with psi_2 instead of psi_1 (order of variables is not changed).

Private Members

Types const *types = nullptr¶: A const pointer to the types of the residues.

Chains const *seqs = nullptr¶: A const pointer to the chain data of the model.

Pairs pairs¶: Local copy of the Verlet list with some pairs potentially excluded (see vlUpdateHook for details).

struct QAContact¶

#include <QuasiAdiabatic.hpp>

A struct detailing the present QA contact.

Public Types

enum Status¶

Enum for the status of a QA contact.

Values:

enumerator FORMING¶: The contact is forming.

enumerator BREAKING¶: The contact is breaking.

enumerator REMOVED¶: The contact has been removed (but since we don’t want to actually remove items from the vector, we just mark it this way.)

Public Members

int i1¶: Index of the first residue.

int i2¶: Index of the second residue.

Status status¶: Status of the contact.

Stats::Type type¶: Type of contact (BB, BS, SB, SS).

double t0¶: Time of either formation or when contact has started breaking.

struct QAFreePair¶

#include <QuasiAdiabatic.hpp>

A struct detailing a free pair that can form a QA contact.

Public Types

enum Status¶

Enum for the status of the free pair.

Values:

enumerator FREE¶: The pair is indeed free.

enumerator TAKEN¶: The pair has formed a QA contact, but since we don’t want to outright remove an element from a list, we only mark it as removed.

Public Members

int i1¶: Index of the first residue.

int i2¶: Index of the second residue.

Status status¶: Status of the pair.

struct QADiff¶

#include <QuasiAdiabatic.hpp>

A “diff” with respect to the QA contact. Since the formation of the contact requires there to be enough “stat slots”, and it’s possible that in an asynchronous execution more contacts would be formed while they shouldn’t, i.e. oversaturate the pool of stat slots, we separate the formation into a “virtual formation” which is done asynchronously and which adds a QADiff to a list of diffs, and a synchronous part which adds these diffs one after another, thus preventing an oversaturation.

Public Functions

inline bool operator<(QADiff const &other) const¶: Linear order on QADiff structures, done in order to ensure the determinism of the process by sorting the diffs’ array.

Public Members

QAContact cont¶: QA contact to add.

int oldIdx¶: Index of a free pair that now got taken.

Stat statDiffs[2]¶: Diffs to Stat that result from the formation of a QA contact.

class QuasiAdiabatic : public mdk::NonlocalForce ¶

#include <QuasiAdiabatic.hpp>

The quasi-adiabatic potential.

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind the force to the simulation. In particular, we fetch the types and chains, load sidechain distances etc.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the computation. In particular, we compute forces between currently-existing pairs and check which of the free pairs could form contacts.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

virtual void syncPart(Dynamics &dynamics) override¶

Synchronous part of the computation. We sort the list of potentially-added contacts and try to add them sequentially.

Parameters: dynamics – Dynamics object to add potential energy and forces to; here it should be unused unless we want to account for the forces and potential energy of the added contacts.

virtual void vlUpdateHook() override¶: Action to be performed when the Verlet list is updated. Here we update the lists of pairs in contact and free pairs, preserving the pairs that were in contact in the old list and is present in the new Verlet list.

Private Functions

virtual vl::Spec spec() const override¶

Generate a spec for the Verlet list.

Returns: Generated spec for the Verlet list.

void computeNH()¶: Compute a list of vectors $n_i$ and $h_i$.

bool geometryPhase(vl::PairInfo const &p, QADiff &diff) const¶

Perform a geometry check between two residues during the formation pass.

Parameters

p – Data regarding the pair currently processed.
diff – A QADiff structure to modify if the check is successful.

Returns

true if the geometry check has succeeded, false otherwise.

Private Members

Types const *types = nullptr¶: Types of residues, used when selecting the sidechain-sidechain potential.

Chains const *chains = nullptr¶: Chains data structure, here used mostly to filter the residues that are at the ends of the residues (there $n_i$ and $h_i$ cannot be computed).

LennardJones bb_lj = LennardJones(5.0 * angstrom, 1.0 * eps)¶: Lennard-Jones potential for the backbone-backbone contacts.

LennardJones bs_lj = LennardJones(6.8 * angstrom, 1.0 * eps)¶: Lennard-Jones potential for the backbone-sidechain contacts.

SidechainLJ ss_ljs[AminoAcid::N][AminoAcid::N]¶: An array of sidechain L-J potentials, one for every pair of types of amino acids.

Vectors n¶: An array of vectors $n_i$, as defined in CPC14.pdf.

Vectors h¶: An array of vectors $h_i$, as defined in CPC14.pdf.

double formationMaxDistSq¶: A maximum distance (across all types of contacts) between residues for the formation to occur. Has no “physical” sense, is used chiefly for optimization.

double hr_abs_min = 0.92¶: Minimum value of <h_i, r_j> for the formation of a backbone contact (with backbone part for i).

double hh_abs_min = 0.75¶: Minimum value of <h_i, h_j> for the formation of a backbone-backbone contact to occur.

double nr_max = 0.5¶: Maximum value of <n_i, r_j> for the formation of a sidechain contact (with sidechain part for i).

double formationTolerance = 1.0¶: A scalar factor by which the “usual” minimum formation distance is multiplied.

double formationTime = 10.0 * tau ¶: A timespan during which a formed QA contact “thermalizes”, i.e. reaches the full depth of the L-J potential. The increase of the depth is linear.

double breakingTolerance = 1.0¶: A scalar factor by which the “usual” maximum distance before breaking is multiplied.

double breakingTime = 10.0 * tau ¶: Timespan during which a broken QA contact “dissipates”, i.e. loses its force. The decrease in the depth of the L-J potential is linear.

std::vector<QAContact> oldPairs¶: The list of old pairs, with which the new pairs are swapped. This is done in order to not have to allocate new memory each time a Verlet list is regenerated.

std::vector<QAContact> pairs¶: A list of current active (or destroyed) QA contacts.

std::vector<QAFreePair> freePairs¶: A list of free pairs, i.e. ones which are not in a QA contact but which potentially could form one.

Stats *stats¶: A reference to Stats variable. Note: this is a non-const reference, and in particular it gets potentially modified in the syncPart.

std::mutex qaDiffsMutex¶: A mutex for the access to the qaDiffs list. We have decided not to employ a separate list for every thread since (we estimate) there shouldn’t be many conflicts to the access; nevertheless we wish to prevent any concurrent access were one to occur.

std::vector<QADiff> qaDiffs¶: A list of QA diffs to pass through in the syncPart, gets new diffs added to in the formation pass.

class Tether : public mdk::Force ¶

#include <Tether.hpp>

Harmonic tether forces between consecutive residues in a chain.

Public Functions

explicit Tether(bool fromNative)¶

Construct a Tether object.

Parameters: fromNative – Whether the equilibrium bond distances should be derived from the native structure or set to a default value (3.8 angstrem).

virtual void bind(Simulation &simulation) override¶

Bind the object to the simulation. Additionally we fetch the chain data to determine which of the pairs (i, i+1) are connected.

Parameters: simulation – Simulation to bind to.

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Private Members

Harmonic harm¶: The underlying harmonic force kernel. It is separated from this class because it’s used in some other places, notably for standard SSBOND interactions.

bool fromNative¶: Whether the equilibrium bond distances should be derived from the native structure or set to a default value (3.8 angstrem).

int n¶: Number of residues, saved for convenience.

Scalars dist0¶: A list of equilibrium distances.

Bytes isConnected¶: isConnected[i] = 1 if pair (i, i+1) is entirely in a chain, i.e. is tethered; 0 otherwise.

class SolidWall : public mdk::Force ¶

#include <SolidWall.hpp>

Solid walls.

Public Functions

virtual void asyncPart(Dynamics &dynamics) override¶

Asynchronous part of the force computation.

Parameters: dynamics – Dynamics object to add potential energy and forces to.

Public Members

Eigen::Hyperplane<double, 3> wall¶

class ExportPDB : public mdk::Hook, private mdk::SimulVar ¶

#include <ExportPDB.hpp>

The hook for exporting current positions of the residues to a PDB file. The positions are stored as MODELs within the PDB file; (at this moment) the file is output only after the resolution of the simulation, i.e. in the destructor. Also, at the moment only the positions (and not for example REMARKs with QA stats or contacts) are saved.

Public Functions

inline explicit ExportPDB(std::filesystem::path modelPath, double period = 1000.0 * microsecond)¶

Construct a ExportPDB hook with a model path and the period.

Parameters

modelPath – Path to the file to which to output the positions.
period – Span between consecutive exports.

~ExportPDB()¶: This destructor saves the to the file.

virtual void bind(Simulation &simulation) override¶

Bind the hook to the simulation. Also, fetch the state.

Parameters: simulation – Simulation to bind to.

virtual void execute(int step_nr) override¶

Execute the hook, i.e. check if the span since the last time the positions were exported is greater than period and if it is, export the positions and set the last time to current time.

Parameters: step_nr – Number of the step of the simulation when the hook is invoked.

Private Members

Model base¶: Reference model to export the current positions to.

pdb::Data data¶: A PDB data object containing the currently stored records (in particular the positions).

int modelIdx = 1¶: Current model serial number (1-indexed).

State const *state = nullptr¶: Const pointer to the state.

std::filesystem::path modelPath¶: Path to the PDB file.

double period¶: Span betwen consecutive exports (in internal time).

double tprev = std::numeric_limits<double>::lowest()¶: Time point when the record was last exported (in internal time). This one is initialized to minimum value so that the first execution always exports.

class Hook¶

#include <Hook.hpp>

This interface defines actions to be performed after the force evaluation, i.e. a hook.

Subclassed by mdk::ExportPDB, mdk::PositionDiff, mdk::ProgressBar

Public Functions

virtual void execute(int step_nr) = 0¶

Execute the hook.

Parameters: step_nr – Number of the step of the simulation when the hook is activated.

class PositionDiff : public mdk::Hook, private mdk::SimulVar ¶

#include <PositionDiff.hpp>

A hook which saves position diffs. I assume this takes the path to the positions as outputted by the Fortran implementation and compares the positions in the consecutive steps with the current positions, whereafter it outputs the diffs to the outputPath.

Public Functions

inline PositionDiff(std::string inputPath, std::string const &outputPath)¶

Create a PositionDiff object with given input and output paths.

Parameters

inputPath – Path to the file with Fortran positions.
outputPath – Path to the diff file to create.

virtual void bind(Simulation &simulation) override¶

Bind the hook to the simulation; also fetch the state and load the reference positions.

Parameters: simulation –

virtual void execute(int step_nr) override¶

Execute the hook, i.e. output the diffs for the current step to the diff file.

Parameters: step_nr – Number of the step of the simulation at which the hook was invoked.

Private Members

std::string inputPath¶: A path to the Fortran reference positions file.

std::ofstream output¶: Position diffs output stream.

std::map<int, Vectors> refPositions¶: Reference positions loaded from the file at inputPath, in the form of a map from the step number to a list of positions.

State const *state = nullptr¶: Saved const pointer to the state of the simulation.

class ProgressBar : public mdk::Hook, private mdk::SimulVar ¶

#include <ProgressBar.hpp>

A progress bar. Chiefly for convenience. Prints a status bar, current (internal and clock) time and total simulation span, and the current potential energy. Note: Because it flushes the output, care must be taken not to invoke it too often.

Public Functions

explicit ProgressBar(double totalTime, double updatePeriod, int width = 70)¶

Construct a ProgressBar object.

Parameters

totalTime – Total span of the simulation; used for normalizing the current time to display it as a percentage.
updatePeriod – How often to update the progrss bar.
width – Width of the progress bar in characters.

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

virtual void execute(int step_nr) override¶

Execute the hook.

Parameters: step_nr – Number of the step of the simulation when the hook is activated.

Private Types

using time_point = std::chrono::high_resolution_clock::time_point¶

Private Members

Simulation *simul = nullptr¶

State *state = nullptr¶

time_point realTime0¶: Real (clock) time when the simulation was started.

double totalTime¶

double updatePeriod¶

double prevTime¶

int width¶

class Harmonic¶

#include <Harmonic.hpp>

Harmonic force “kernel”, i.e. a separated class responsible only for the computation of the formula. The definitions are inlined in order for the compiler to inline them.

Public Functions

Harmonic() = default¶

inline Harmonic(double H1, double H2)¶

inline void computeV(double dx, double &V, double &dV_dx) const¶

Compute the potential energy of the force field.

Parameters

dx – Displacement, i.e. the difference between current length and the equilibrium length
V – Variable to add the potential to.
dV_dx – Variable to add the derivative to.

template<typename T1, typename T2> inline void computeF(VRef unit, double dx, double &V, T1 F1, T2 F2) const¶

Compute and add the harmonic force between two residues. The templates are here in order for us to be able to pass Eigen expressions to it.

Template Parameters

T1 – Type of an lvalue to add the force on the first residue to.
T2 – Type of an lvalue to add the force on the second residue to.

Parameters

unit – Normalized vector between the residues.
dx – Displacement, i.e. the difference between current length and the equilibrium length;
V – Variable to add the potential to.
F1 – Lvalue to add the force on the first residue to.
F2 – Lvalue to add the force on the second residue to.

Public Members

double H1 = 50.0 * eps / pow(angstrom, 2.0)¶: “Harmonic” part of the harmonic potential.

double H2 = 0.0¶: “Anharmonic” part of the harmonic potential. From what I have seen it’s pretty much always zero, perhaps it’s a legacy/deprecated feature.

class LennardJones¶

#include <LennardJones.hpp>

Standard Lennard-Jones potential.

Public Functions

LennardJones() = default¶

inline LennardJones(double r_min, double depth)¶

inline double cutoff() const¶

The cutoff distance for the L-J potential. I’m not sure what is the “industry standard” cutoff but at 2.5 sigma the value of the potential is -0.016 epsilon, which seems fine as a cutoff. For comparison, at 2 sigma the value is -0.06225 epsilon which seems too high to cut off.

Returns: The cutoff distance.

inline void computeV(double norm, double &V, double &dV_dn) const¶

Compute the potential energy of the force field.

Parameters

norm – Distance between the residues.
V – Variable to add the potential to.
dV_dn – Variable to add the derivative to.

template<typename T1, typename T2> inline void computeF(VRef unit, double norm, double &V, T1 F1, T2 F2) const¶

Compute and add the L-J force between two residues. The templates are here in order for us to be able to pass Eigen expressions to it.

Template Parameters

T1 – Type of an lvalue to add the force on the first residue to.
T2 – Type of an lvalue to add the force on the second residue to.

Parameters

unit – Normalized vector between the residues.
norm – Distance between the residues.
V – Variable to add the potential to.
F1 – Lvalue to add the force on the first residue to.
F2 – Lvalue to add the force on the second residue to.

Public Members

double r_min = 5.0 * angstrom ¶

double depth = 1.0 * eps ¶

class ShiftedTruncatedLJ¶

#include <ShiftedTruncatedLJ.hpp>

A shifted and truncated version of the L-J potential, used chiefly by the Pauli exclusion force field.

Public Functions

ShiftedTruncatedLJ() = default¶

inline ShiftedTruncatedLJ(double r_cut, double depth)¶

inline double cutoff() const¶

inline void computeV(double norm, double &V, double &dV_dn) const¶

Compute the potential energy of the force field.

Parameters

norm – Distance between the residues.
V – Variable to add the potential to.
dV_dn – Variable to add the derivative to.

template<typename T1, typename T2> inline void computeF(VRef unit, double norm, double &V, T1 F1, T2 F2) const¶

Compute and add the L-J force between two residues. The templates are here in order for us to be able to pass Eigen expressions to it.

Template Parameters

T1 – Type of an lvalue to add the force on the first residue to.
T2 – Type of an lvalue to add the force on the second residue to.

Parameters

unit – Normalized vector between the residues.
norm – Distance between the residues.
V – Variable to add the potential to.
F1 – Lvalue to add the force on the first residue to.
F2 – Lvalue to add the force on the second residue to.

Public Members

double r_cut = 4.0 * angstrom ¶

double depth = 1.0 * eps ¶

class SidechainLJ¶

#include <SidechainLJ.hpp>

Sidechain version of the L-J potential (for a single pair of depth and sink_max). In the CPC14.pdf it’s described as an L-J potential with a well bounded from both ends; in practice this is implemented as a one-sided well, i.e. for r < sink_max, V = - depth and the other side of the well is handled by the Pauli exclusion.

Public Functions

SidechainLJ() = default¶

inline SidechainLJ(double depth, double sink_max)¶

inline double cutoff() const¶

inline void computeV(double norm, double &V, double &dV_dn) const¶

Compute the potential energy of the force field.

Parameters

norm – Distance between the residues.
V – Variable to add the potential to.
dV_dn – Variable to add the derivative to.

template<typename T1, typename T2> inline void computeF(VRef unit, double norm, double &V, T1 F1, T2 F2) const¶

Compute and add the L-J force between two residues. The templates are here in order for us to be able to pass Eigen expressions to it.

Template Parameters

T1 – Type of an lvalue to add the force on the first residue to.
T2 – Type of an lvalue to add the force on the second residue to.

Parameters

unit – Normalized vector between the residues.
norm – Distance between the residues.
V – Variable to add the potential to.
F1 – Lvalue to add the force on the first residue to.
F2 – Lvalue to add the force on the second residue to.

Public Members

double depth = 1.0 * eps ¶

double sink_max = 10.0 * angstrom ¶

class Model¶

#include <Model.hpp>

A model object, which serves as one of two objects (other being param::Parameters) which instantiate the simulation. Contains residues (with types, messes, positions etc.), chains and structured parts, topological information (i.e. the box shape), along with an array of utilities for the modification of the model, in particular initializing the velocities or preparing the model into a conformation (a line or a self-avoiding walk).

Public Functions

Residue &addResidue(Chain *chain = nullptr)¶

Chain &addChain()¶

Contact &addContact()¶

StructuredPart &addSP()¶

void morphIntoLine()¶

void morphIntoSAW(Random &rand, bool useTop = false, double density = 1e-4 * atom / pow(angstrom, 3.0), double minDist = 4.56 * angstrom)¶

Morphing the model into a self-avoiding walk. A custom version with clarity in mind, doesn’t replicate the positions from the Fortran code which makes testing difficult.

Note: it doesn’t use native bond angles correctly at the moment.

Parameters

rand – Random number source to use.
useTop – Whether to create and use the topology during the creation process.
density – Target density of the box. 0 means an infinite box?
minDist – Minimum distance for the conformation to be considered intersecting, and therefore invalid.

void legacyMorphIntoSAW(Random &rand, bool useTop = false, double density = 1e-4 * atom / pow(angstrom, 3.0), double minDist = 4.56 * angstrom, bool nativeBondLen = false)¶

Morphing the model into a self-avoiding walk. A version which replicates the behavior of the Fortran code but in a more clear way.

Parameters

rand – Random number source to use.
useTop – Whether to create and use the topology during the creation process.
density – Target density of the box. 0 means an infinite box?
minDist – Minimum distance for the conformation to be considered intersecting, and therefore invalid.
nativeBondLen – Whether to derive the bond lengths from the native structure. By default the bond lengths are 3.8 angstrem; if nativeBondLen is true, the bond length is taken to be a mean bond length in the native structure.

void exactLegacySAW(Random &rand, bool useTop = false, double density = 0.0 * atom / pow(angstrom, 3.0), double minDist = 4.56 * angstrom, double bond = 3.8 * angstrom)¶

Morphing the model into a self-avoiding walk. A version which replicates (one-by-one) the code of the Fortran version. It is also less intuitive because of it.

Parameters

rand – Random number source to use.
useTop – Whether to create and use the topology during the creation process.
density – Target density of the box. 0 means an infinite box?
minDist – Minimum distance for the conformation to be considered intersecting, and therefore invalid. Note, from what I see the intersection is checked only for each chain internally, and not across the chains.
bond – Default bond length to use.

void initVelocity(Random &rand, double temperature, bool useMass = true)¶

Initialize velocities of the residues, assuming they are placed in a medium with a given temperature.

Parameters

rand – Random number source to use.
temperature – Temperature of the medium.
useMass – Whether to account for different masses of the residues.

StructuredPart &addContactMap(cmap::ContactMap const &contactMap)¶

void addCMapContacts(cmap::ContactMap const &contactMap, Chain &chain)¶

void useAverageMasses()¶: Set the masses of all residues to be equal to the mean of the current masses.

Public Members

int n = 0¶: The number of residues in the model.

std::vector<Residue> residues¶

std::vector<Chain> chains¶

std::vector<Contact> contacts¶

std::vector<StructuredPart> structuredParts¶

Topology top¶: Topology, or the box size.

Private Functions

std::vector<std::pair<Residue*, Residue*>> nonlocalPairs()¶

Generates the pairs of residues separated by at least two bonds, used in the SAW procedures.

Returns: Pairs of residues.

struct Chain¶

Public Members

int idx¶

int start¶

int end¶

std::vector<int> structuredParts¶: A list of indices of structured parts that are overlayed over the chain.

struct Contact¶

Public Members

int idx¶

int res[2]¶

ContactType type¶

double dist0¶

struct Residue¶

#include <Model.hpp>

A structure representing a single “residue”, or a pseudoatom, of the model.

Public Members

int idx¶: Index of the residue in residues.

int chainIdx¶: Index of the parent chain, or -1 if it doesn’t belong to any chain.

Vector r¶

Vector v¶

double mass¶

ResType type¶

std::optional<Vector> nat_r¶: Native positions of the residues.

struct StructuredPart¶

Public Members

int idx¶

int off¶

int len¶

std::vector<double> angle¶

std::vector<double> dihedral¶

class Simulation¶

#include <Simulation.hpp>

The main class of the library, responsible for:

storing the simulation state;
storing the forces and other associated objects, and (optionally) binding them to the simulation;
running, or rather stepping through, the simulation.

Public Functions

inline Simulation(Model model, param::Parameters params)¶

Initialize simulation from a Model object and the parameters.

Parameters

model – Model of the simulation.
params – Parameters of the simulation.

template<typename Data> inline Data const &data()¶

Access the “static” data as created by the DataFactory. We access the data by the type, but insofar as these objects represent functions of the state and the simulation, this is (in my estimation) appropriate. Providing an invalid type will lead to linker error.

Template Parameters: Data – Type of static data to access.
Returns: Const reference to the data object.

template<typename Var> inline Var &var()¶

Accesses the variables associated with the simulation. The access is by type, thus it’s most suited to variables that are effectively functions of the simulation (for example the physical state or the Verlet list). If the referenced variable has not been yet added, the function tries to emplace it (if it’s default-constructible) such is, for example, the case with State or vl::List. Cannot be run after simulation initialization. The reference is never invalidated.

Template Parameters: Var – Type of the variable to access.
Returns: Accessed variable. The function returns a non-const reference, thus care must be taken as to modify the variables only on agreed-upon points of the simulation (say, integration of forces).

template<typename T, typename ...Args> inline T &add(Args&&... args)¶

Add a variable to the simulation. Depending on the type of the variable, additional actions may be performed (such as binding to the simulation or adding forces and/or nonlocal forces to appropriate lists). Cannot be run after simulation initialization. The reference is never invalidated.

Template Parameters

T – Type of the variable to add.
Args – Types of arguments to use in construction.

Parameters

args – Values of arguments to use in construction.

Returns

Reference to a constructed and emplaced T(Args…).

inline void addAsyncTask(std::function<void()> const &f)¶

Add an “async task”. Cannot be run after simulation initialization.

Parameters: f – Function/lambda to add. (Seemingly the lambda cannot capture the environment, i.e. must be trivial)

void init()¶: Initialize the simulation.

void step()¶: Step once through the simulation; the time step is as determined by the used integrator (usually 5 ps).

void step(double t)¶

Step a total of time t through the simulation.

Parameters: t – Time to advance.

Private Functions

void calcForces()¶: Internal function for invoking force fields.

Private Members

Model model¶: Model of the simulation.

param::Parameters params¶: Parameters of the simulation.

DataFactory df¶: Data factory.

std::unordered_map<std::type_index, std::shared_ptr<void>> vars¶: A heterogeneous collection of objects (stored as shared pointers to void) indexed by their types.

std::vector<Force*> forces¶: Force fields of the simulation.

std::vector<NonlocalForce*> nonlocalForces¶: List of nonlocal forces of the simulation - the distinction is used when the Verlet list is updated, so as to invoke appropriate update functions.

std::vector<Hook*> hooks¶: Hooks of the simulation.

std::vector<std::function<void()>> asyncTasks¶: Async tasks.

Integrator *integrator = nullptr¶: Pointer to the integrator used in the simulation.

bool initialized = false¶: Whether simulation has yet been initialized; used for running the init function only once.

int step_nr = 0¶: Number of steps performed.

class SimulVar¶

#include <SimulVar.hpp>

A “simulation variable”, i.e. an object which can be bound to a particular simulation (usually to finish instantiation).

Subclassed by mdk::ExportPDB, mdk::Force, mdk::HeuresticBA, mdk::HeuresticDihedral, mdk::Integrator, mdk::NativeBA, mdk::NativeDihedralBase, mdk::PositionDiff, mdk::ProgressBar, mdk::State, mdk::Stats, mdk::vl::List

Public Functions

virtual void bind(Simulation &simulation) = 0¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

struct Stat¶

#include <Stat.hpp>

An object containing the coordination stats for a given residue. As opposed to the Fortran version, where the numbers represent the number of contacts formed, here it’s the number of slots remaining. This simplifies the checking procedure, as we can simply compare the coordination numbers with zero. We also provide facilities for manipulating these stats, such as addition etc. This object can also naturally represent a stat difference.

Public Functions

inline bool valid() const¶

Checks whether a stat is valid; in particular, one can check whether the formation of a contact will exceed the number of available slots this way.

Returns: Whether a stat is valid.

inline Stat &operator+=(Stat const &other)¶

inline Stat operator+(Stat const &other) const¶

inline Stat operator-() const¶

inline Stat &operator-=(Stat const &other)¶

inline Stat operator-(Stat const &other) const¶

Public Members

int8_t backbone = 0¶: Number of remaining backbone contacts that can be formed.

int8_t sidechain = 0¶: Number of remaining sidechain contacts that can be formed.

int8_t hydrophobicSS = 0¶: Number of remaining sidechain contacts that can be formed with hydrophobic residues.

int8_t polarSS = 0¶: Number of remaining sidechain contacts that can be formed with polar residues.

class Stats : public mdk::SimulVar ¶

#include <Stats.hpp>

An object containing stats for the residues. Also initializes these numbers in one place.

Public Types

enum Type¶

Type of a contact (backbone-backbone etc.).

Values:

enumerator BB¶

enumerator BS¶

enumerator SB¶

enumerator SS¶

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

inline void creationDiffs(int i1, int i2, Type type, Stat diffs[2]) const¶

Generate diffs for if a contact between residues i1 and i2 of type type were created. It’s inline so as for the compiler to insert it and not have to call the function.

Parameters

i1 – First residue.
i2 – Second residue.
type – Type of a contact.
diffs – An array of diffs; diffs[0] is the diff for the first residue, diffs[1] is the diff for the second residue.

Public Members

std::vector<Stat> stats¶: List of stats for each residue.

Private Members

Types const *types = nullptr¶: Types of amino acids.

std::vector<param::Polarization> polarization¶: List of polarization values taken from the param::Parameters object.

class Integrator : public mdk::SimulVar ¶

#include <Integrator.hpp>

An integrator object; after the forces have been computed, it takes the forces and adjusts the positions, velocities etc.

Subclassed by mdk::LangPredictorCorrector, mdk::Leapfrog

Public Functions

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

virtual void init() = 0¶: Initialize the state of the integrator.

virtual void integrate() = 0¶

Protected Attributes

State *state = nullptr¶

class LangPredictorCorrector : public mdk::Integrator ¶

#include <LangPredictorCorrector.hpp>

Langevin predictor-corrector fifth-order integrator, combined with a Langevin noise. Apparently these two are combined because in the Fortran code, the velocities are set in lang or lang_mass. Now, this looks like a bug, but we replicate it nonetheless.

Public Functions

inline explicit LangPredictorCorrector(double dt)¶

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

virtual void init() override¶: Initialize the state of the integrator.

virtual void integrate() override¶

Public Members

double gamma = 2.0 * f77mass / tau ¶: Value of gamma for the Langevin noise.

Private Functions

void generateNoise()¶

Private Members

double dt¶

Masses m¶

Vectors y0¶

Vectors y1¶

Vectors y2¶

Vectors y3¶

Vectors y4¶

Vectors y5¶

Random *random = nullptr¶

std::vector<Random> rngs¶

Vectors gaussianNoise¶

bool initialized = false¶

double temperature = 0.35 * eps_kB ¶

class Leapfrog : public mdk::Integrator ¶

#include <Leapfrog.hpp>

Leapfrog integration algorithm.

Public Functions

inline explicit Leapfrog(double dt)¶

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

virtual void init() override¶: Initialize the state of the integrator.

virtual void integrate() override¶

Private Members

double dt = 0.005 * tau ¶

Masses m¶

Vectors a_prev¶

struct Dynamics¶

#include <State.hpp>

An object containing the dynamical state of the simulation, i.e. forces and the potential energy.

Public Functions

inline void zero(int n)¶

Public Members

double V = 0.0¶

Vectors F¶

class State : public mdk::SimulVar ¶

#include <State.hpp>

An object containing the physical state of the simulation, i.e. current positions, velocities of the residues, time and the box shape.

Public Functions

void prepareDyn()¶

void updateWithDyn(Dynamics const &othDyn)¶

virtual void bind(Simulation &simul) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

void exportTo(Model &model) const¶

Set the positions and velocities of the residues in a Model object to the current values.

Parameters: model –

Public Members

int n¶

Vectors r¶

Vectors v¶

double t¶

Topology top¶

Dynamics dyn¶

struct AAAtomInfo¶

#include <AminoAcid.hpp>

Info pertaining to a single atom type in an amino acid (for example SG in cysteine).

Public Members

double radius¶

bool inBackbone = false¶

struct AminoAcidInfo¶

#include <AminoAcid.hpp>

Info pertaining to an amino acid type.

Public Members

double mass¶

std::unordered_set<std::string> heavyAtoms¶: A set of names of the heavy atoms for the amino acid.

std::unordered_map<std::string, AAAtomInfo> atomInfo¶: A map from the heavy atom names to the atom info structures.

class AminoAcid¶

#include <AminoAcid.hpp>

An object representing an amino acid type, along with a number of facilities, conversion from/to other types, general static data etc.

Public Functions

AminoAcid() = default¶

inline explicit constexpr AminoAcid(AminoAcidIdx type)¶

explicit operator AminoAcidIdx const&() const¶

inline constexpr AminoAcid(int8_t x)¶

operator int8_t() const¶

explicit AminoAcid(char code)¶

explicit operator char const&() const¶

explicit AminoAcid(std::string const &name)¶

explicit operator std::string() const¶

bool operator==(AminoAcid const &other) const¶

bool operator!=(AminoAcid const &other) const¶

bool operator<(AminoAcid const &other) const¶

bool operator<=(AminoAcid const &other) const¶

bool operator>(AminoAcid const &other) const¶

bool operator>=(AminoAcid const &other) const¶

AminoAcidInfo const &info() const¶

Public Static Functions

static std::vector<AminoAcid> aminoAcids()¶

static std::vector<AminoAcidIdx> types()¶

static std::string codes()¶

static bool isProper(std::string const &name)¶

Check whether a string corresponds to an amino acid proper name.

Parameters: name – Text to check.
Returns: true if name corresponds to an amino acid, false otherwise.

static std::vector<std::string> names()¶

Public Static Attributes

static constexpr const int N = 20¶

Private Members

AminoAcidIdx type¶

Friends

friend struct std::hash< AminoAcid >

class ContactType¶

#include <ContactType.hpp>

An object representing a contact type, along with a number of facilities, conversions from/to other types etc.

Public Functions

ContactType() = default¶

inline explicit constexpr ContactType(ContactTypeIdx code)¶

explicit operator ContactTypeIdx const&() const¶

inline constexpr ContactType(int8_t x)¶

operator int8_t() const¶

explicit ContactType(std::string const &name)¶

explicit operator std::string() const¶

Private Members

ContactTypeIdx code¶

Friends

friend struct std::hash< ContactType >

class Random¶

#include <Random.hpp>

A random number generator. We use two versions: a legacy version taken from Fortran, and a modern version. Aside from sampling from [0, 1], we add sampling from [a, b], N(0, 1), N(mu, sigma^2) or points on S^2.

Public Functions

inline Random(int seed)¶

inline double uniform()¶

inline Random getNewRandom()¶

Random(Random const &oth) = default¶

inline double uniform(double a, double b)¶

inline double normal()¶

inline std::pair<double, double> two_normals()¶

inline double normal(double mu, double sigma)¶

inline Eigen::Vector3d sphere()¶

Private Functions

Random() = default¶

Private Members

uint64_t state¶

class ResType¶

#include <ResType.hpp>

An object representing a residue type, along with a number of facilities, conversions from/to other types etc. At this moment, it is equivalent with the amino acids, but one may want to add other pseudo-atoms, sites, water molecules or groups etc. so we thought it wiser to separate the two.

Public Functions

ResType() = default¶

inline explicit constexpr ResType(ResTypeIdx code)¶

explicit operator ResTypeIdx const&() const¶

inline constexpr ResType(int8_t x)¶

operator int8_t() const¶

explicit ResType(std::string const &name)¶

explicit operator std::string() const¶

double mass() const¶

Private Members

ResTypeIdx code¶

Friends

friend struct std::hash< ResType >

class Topology¶

#include <Topology.hpp>

An object representing the “topology” of the simulation box. In particular, we want to compute PBC-adjusted distances between residues.

Public Functions

void setCell(VRef cell)¶

inline void fix(Vector &v) const¶

inline Vector operator()(Vector v) const¶

Public Members

Vector cell¶

Vector cellInv¶

bool use[3] = {false, false, false}¶

namespace cmap¶

class ContactMap¶

#include <ContactMap.hpp>

A POD representing a contact map (a.k.a. structured part), as represented in the legacy contact map file): specifically, the location and extent of the structured part, associated native bond and dihedral angles and Go-model contacts.

This data structure doesn’t provide its own means of parsing (cmap::LegacyParser must be used, for example), since a change of format shouldn’t influence this POD.

Public Members

int len¶: Length of the structured part.

int offset¶: Offset of the structured part (w.r.t. the chain beginning).

std::vector<Contact> contacts¶: List of Go-model contacts.

std::vector<double> angle¶: angle[i] denotes the native bond angle value for triple (i-1, i, i+1), relative to the position as indicated by /p off, i.e. if (i+off-1, i+off, i+off+1) relative to the start of the chain. Note: the native angles near the ends of the structured part are undefined.

std::vector<double> dihedral¶: dihedral[i] denotes the native dihedral angle value for quadruple (i-2, i-1, i, i+1), relative to the position as indicated by /p off, i.e. if (i+off-2, i+off-1, i+off, i+off+1) relative to the start of the chain. Note: the native angles near the ends of the structured part are undefined.

struct Contact¶

Public Members

int res[2]¶: Residues forming the contact.

double dist0¶: Native length of the contact.

class LegacyParser¶

#include <LegacyParser.hpp>

A parser of contact map files in the legacy format (as described in legacy README.txt).

Public Functions

ContactMap read(std::istream &is)¶

Read contact map from a stream.

Parameters: is – Input stream wherefrom to read the contact map
Returns: Parsed contact map.

std::ostream &write(std::ostream &os, ContactMap const &cmap)¶

Write contact map to a stream.

Parameters

os – Output stream whereto to write the contact map
cmap – Contact map to write

Returns

os

namespace param¶

Enums

enum Polarization¶

Polarization values, as derived from CPC14.pdf and from reading cg.f.

Values:

enumerator GLY_PRO¶

enumerator HYDROPHOBIC¶

enumerator MISSING¶

enumerator POLAR¶

enumerator POLAR_NEG¶

enumerator POLAR_POS¶

class LegacyParser¶

#include <LegacyParser.hpp>

A parser of parameter files in the legacy format (as described in legacy README.txt).

Public Functions

Parameters read(std::istream &is)¶

Read parameter file from a stream.

Parameters: is – Input stream wherefrom to read the contact map
Returns: Parsed parameter file.

std::ostream &write(std::ostream &os, Parameters const &data)¶

Write parameter file to a stream.

Parameters

os – Output stream whereto to write the contact map
data – Parameter file to write

Returns

os

class Parameters¶

#include <Parameters.hpp>

A POD representing the parameters from the parameter file. Arguably the legacy file could be split so that everything wouldn’t be in one structure.

Public Types

using Coeffs = std::vector<double>¶

template<typename T> using PerPTData = std::unordered_map<PairType, T>¶

template<typename T> using PerAcidData = std::unordered_map<AminoAcid, T>¶

template<typename T> using PerPairData = std::unordered_map<std::pair<AminoAcid, AminoAcid>, T>¶

Public Functions

Parameters()¶

Public Members

Coeffs defAngleParams¶: “Default” angle parameters. These are I assume for the simplified angle potential (option lsimpang).

PerPTData<Coeffs> angleParams¶: A set of bond angle parameters (for heurestic potential) for each pair of types of residues i, i+1 in a triple (i-1, i, i+1).

PerPTData<Coeffs> dihedralParams¶: A set of dihedral angle parameters (for heurestic potential) for each pair of types of residues i-1, i in a quadruple (i-2, i-1, i, i+1).

PerAcidData<SpecificityParams> specificity¶: Specificity parameters for every amino acid type.

PerAcidData<double> radius¶: Radii of amino acids (for use in non-full-atomic derivation of a contact map from the native structure).

PerPairData<double> pairwiseMinDist¶: A maximum in a distribution of sidechain-sidechain contacts between two given amino acids.

std::optional<PerPairData<double>> mjMatrix¶: An optional MJ matrix of Lennard-Jones potential depths depending on the amino acid types in the pair.

struct SpecificityParams¶

#include <Parameters.hpp>

Specificity parameters POD.

Public Members

Polarization polarization¶: Polarization of an amino acid.

int maxSidechain¶: Maximum number of sidechain contacts an amino acid can form.

int maxHydrophobicSS¶: Maximum number of sidechain contacts an amino acid can form with hydrophobic amino acids.

int maxPolarSS¶: Maximum number of sidechain contacts an amino acid can form with polar amino acids.

namespace pdb¶

Enums

enum Mode¶

A mode enum, which modifies how to position and parse strings. Note: most enums in the program are enum classes; here we need a standard enum in order for combining via | to work

Values:

enumerator Exact¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

enumerator Trim¶: When parsing, the final value is a trimmed string read from the line.

enumerator Left¶: When writing, the string is written left-aligned in the line.

enumerator Right¶: When writing, the string is written right-aligned in the line.

Functions

Data &operator<<(Data &data, records::Record const &record)¶

Add a record to the PDB data object.

Parameters

data – PDB data object.
record – PDB record to add.

Returns

data but with added record

Data &operator<<(Data &data, Data const &other)¶

Append a list of records from another PDB data object.

Parameters

data – PDB data object
other – Another PDB data object from where to get the records

Returns

data but with added records from other

class Data¶

#include <Data.hpp>

PDB file in the “raw format”, i.e. as a list of PDB records. Certain utility functions, like parsing from/to pdb::Model and retrieving only atom positions, are provided.

Public Functions

Data() = default¶

explicit Data(pdb::Model const &model)¶

Convert a PDB model into a set of records. In particular, extracts atom positions into ATOM records, contacts into SSBOND and LINK records, cell shape into CRYST1.

Parameters: model – PDB model to convert.

pdb::Model asModel() const¶

Convert the data into a PDB model (ATOM records into atoms, SSBOND and LINK records into contacts, CRYST1 into cell shape).

Returns: PDB model generated from the data.

Public Members

std::vector<records::Record> records¶: List of PDB records in the PDB file.

Public Static Functions

static Data onlyAtoms(pdb::Model const &model)¶

Extract only ATOM records from a PDB model. Useful for outputting a single file with lots of conformations of the same amino acid chain, for example when we output PDB model with consecutive positions in the simulation.

Parameters: model – PDB model to convert.
Returns: PDB data object containing only ATOM records.

Friends

friend Data &operator<<(Data &data, records::Record const &record)¶

Add a record to the PDB data object.

Parameters

data – PDB data object.
record – PDB record to add.

Returns

data but with added record

friend Data &operator<<(Data &data, Data const &other)¶

Append a list of records from another PDB data object.

Parameters

data – PDB data object
other – Another PDB data object from where to get the records

Returns

data but with added records from other

class Model¶

#include <Model.hpp>

PDB model object. As opposed to Model, here we store atomic data, chiefly for the use in the derivation of full-atomic contact maps. Also, non-full-atomic contact maps may be derived; Model doesn’t contain such a method because it’s ill-defined for models obtained from Sequence file, which doesn’t have native residue positions.

Public Functions

Atom &addAtom(int serial, Residue *res = nullptr)¶

Add an atom to a residue.

Parameters

serial – Serial number of an atom
res – Parent residue, or nullptr

Returns

Reference to the added atom.

Residue &addResidue(int serial, Chain *chain = nullptr)¶

Add a residue to a chain.

Parameters

serial – Serial number of a residue to add.
chain – Parent chain or nullptr

Returns

Reference to the added residue.

Chain &addChain(char serial)¶

Add a chain with a specified serial character.

Parameters: serial – Serial character of the chain.
Returns: Reference to the added chain.

Contact &addContact()¶

Add a new contact.

Returns: Reference to the added contact.

Model() = default¶

explicit Model(mdk::Model const &coarse)¶

Construct a pdb::Model from a Model. Doesn’t extrapolate extra atoms, i.e. only CA atoms are in the final model.

Parameters: coarse – Model from which to construct the PDB model.

void addContactsFromAtomOverlap()¶: Adds contacts based on the overlap of atoms. This method also derives types of contacts based on whether the atoms are in the backbone or the sidechain of an amino acid.

void addContactsFromResOverlap(param::Parameters const &params)¶

Adds contacts based on approximate residue overlap, given the aminoAcidRadii data from the parameter file.

Parameters: params – Parameter file with amino acid radii data.

mdk::Model coarsen()¶

Construct a mdk::Model from the PDB model. Collapses residues to a single Model residue and copies the contacts (if two atoms are connected, the residues containing them are connected in the final model).

Returns: A reduced (“coarsened”) model.

Public Members

std::unordered_map<int, Atom> atoms¶: A map from the serial number of an atom to the Atom structure. Note that the order of serial numbers is (1) not necessarily consecutive, (2) not starting from 0, in some degenerate cases. Also, adding new atoms doesn’t invalidate pointers or references.

std::unordered_map<int, Residue> residues¶: A map from residue serial number to Residue struct. Note that the order of serial numbers is (1) not necessarily consecutive, (2) not starting from 0, in some degenerate cases. Also, adding new residues doesn’t invalidate pointers or references.

std::unordered_map<char, Chain> chains¶: A map from the chain serial character to the chain.

std::vector<Contact> contacts¶: A vector of contacts. Note that references are (potentially) invalidated when adding new contacts.

Topology top¶: “Topology”, i.e. box shape for PBC.

struct Atom¶

#include <Model.hpp>

A struct containing data for a single atom.

Public Members

int serial¶: PDB serial number of an atom. Zero-indexed.

int idxInRes¶: Index of an atom in the Residue::atoms list of the parent residue, or -1 if it’s not a part of any residue.

Residue *res¶: Parent residue.

Eigen::Vector3d r¶: (Native) position of the atom.

std::string type¶: Type of an atom (in the std::string form).

struct Chain¶

#include <Model.hpp>

A struct containing PDB chain data.

Public Members

char serial¶: Serial (alphanumeric) character of a chain.

std::vector<Residue*> residues¶: Constituent resides of the chain. Ordered.

struct Contact¶

#include <Model.hpp>

A struct containing contacts. Technically this isn’t represented 1-to-1 in the PDB data, but an unification of SSBONDs and LINKs is useful, especially for the reduction of the pdb::Model into Model.

Public Members

int idx¶: Index of the contact in contacts vector.

Atom *atom[2]¶: Pointers to atoms that comprise the contact.

std::string type¶: Type of contact. Takes (at this moment) values “SSBOND”, “NAT”, “NAT_BB”, “NAT_BS”, “NAT_SB”, “NAT_SS”.

double dist0¶: Native contact length.

struct Residue¶

#include <Model.hpp>

A struct containing data for a PDB residue.

Public Functions

Atom *find(std::string const &type)¶

Find an atom belonging to the residue by the type. (Used primarily for finding CG atoms in SSBONDs).

Parameters: type – Atom type to find.
Returns: Pointer to found atom, or nullptr if such an atom has not been found.

Public Members

int serial¶: Serial number of a residue. 0-indexed.

int idxInChain¶: Index in the Chain::residues of the parent chain, or -1 if it’s not in any chain.

Chain *chain¶: Parent chain.

std::string type¶: Type of residue (in std::string form)

std::vector<Atom*> atoms¶: List of constituent atoms for the residue.

class Parser¶

#include <Parser.hpp>

A parser of PDB files.

Public Functions

Parser()¶

Data read(std::istream &is)¶

Read PDB data from an input stream.

Parameters: is – Input stream wherefrom to read the data.
Returns: Read PDB data.

std::ostream &write(std::ostream &os, Data const &data)¶

Write formatted PDB data to an output stream.

Parameters

os – Output stream whereto to write the PDB data.
data – PDB data to write.

Returns

os.

Private Members

std::unordered_map<int, std::shared_ptr<RecordParser>> parsers¶: (Internal) set of record parsers.

class Field¶

#include <Field.hpp>

An object representing a single field of a PDB record, along with the means of parsing and printing it. The field objects signal parsing failure by throwing.

Subclassed by mdk::pdb::Char, mdk::pdb::Integer, mdk::pdb::Literal, mdk::pdb::Real, mdk::pdb::String, mdk::pdb::SymOp

Public Functions

virtual void read(const std::string &s) = 0¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const = 0¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

class Integer : public mdk::pdb::Field ¶

#include <Field.hpp>

Integer field.

Public Functions

inline Integer(int i, int j, int &v, int offset = 0)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int j¶: End of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int offset¶: A value that is added to the read value and subtracted when a value is written. It’s chiefly used to convert 0-indexed ordinals to 1-indexed and vice versa.

int *v¶: A pointer to the integer to save the value to and load the value from. It allows us in particular to have the same field class for many different PDB records, only having to provide different pointer to the underlying value.

class Real : public mdk::pdb::Field ¶

#include <Field.hpp>

Real (floating number) field.

Public Functions

inline Real(int i, int j, int n, int m, double &v, double scalar = 1.0)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int j¶: End of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int n¶: Width of the read real number. Follows the conventions from the PDB format docs and Fortran specifiers (it’s n in f{n}.{m})

int m¶: Precision of the read real number. Follows the conventions from the PDB format docs and Fortran specifiers (it’s m in f{n}.{m})

double *v¶: A pointer to the double to save the value to and load the value from. It allows us in particular to have the same field class for many different PDB records, only having to provide different pointer to the underlying value.

double scalar¶: Scalar by which to multiply a read value and divide a written value. In particular, it allows us to convert PDB file formats to internal units and back.

class String : public mdk::pdb::Field ¶

#include <Field.hpp>

String field. By default read values are trimmed, and written values are right-aligned.

Public Functions

inline String(int i, int j, std::string &v, int8_t mode = Trim | Right)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int j¶: End of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

std::string *v¶: A pointer to the string to save the value to and load the value from. It allows us in particular to have the same field class for many different PDB records, only having to provide different pointer to the underlying value.

int8_t mode¶: Mode of string parsing and writing. See Mode for details.

class Char : public mdk::pdb::Field ¶

#include <Field.hpp>

(Single) character field.

Public Functions

inline Char(int i, char &v)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: The location of the character in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

char *v¶: A pointer to the character to save the value to and load the value from. It allows us in particular to have the same field class for many different PDB records, only having to provide different pointer to the underlying value.

class SymOp : public mdk::pdb::Field ¶

#include <Field.hpp>

“Symmetry operator” field. At the moment it’s parsed just like a normal string.

Public Functions

inline SymOp(int i, int j, std::string &v)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int j¶: End of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

std::string *v¶: A pointer to the string to save the value to and load the value from. It allows us in particular to have the same field class for many different PDB records, only having to provide different pointer to the underlying value.

class Literal : public mdk::pdb::Field ¶

#include <Field.hpp>

A literal. As opposed to a normal field, here there’s no reading or writing values by pointers, rather comparing the read value with the literal and printing out the literal. It’s chiefly used for parsing the record headers, like “ATOM ” or “MODEL “.

Public Functions

inline Literal(int i, int j, std::string lit)¶

virtual void read(const std::string &s) override¶

Read the value from a string (usually a line of a PDB file).

Parameters: s – String to read the value from.

virtual void write(std::string &s) const override¶

Write the value to a string (usually a line of a PDB file). Overflows are handled depending on the type of field.

Parameters: s – String to write the value to.

Private Members

int i¶: Beginning of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

int j¶: End of the field in the line. Inclusive and 1-indexed, as in the PDB format docs, in order to make writing new fields easier.

std::string lit¶: Literal to compare the read value with and print out to the line when writing.

class RecordParser¶

#include <RecordParsers.hpp>

A parser of a single record. Because of the “bidirectionality” of Field structures and ability to swap out the pointers, we can elegantly represent as a list of the fields in the line of a record, and the parsing and writing of the record are essentially provided out-of-the-box.

Subclassed by mdk::pdb::AtomParser, mdk::pdb::Cryst1Parser, mdk::pdb::EndParser, mdk::pdb::EndmdlParser, mdk::pdb::HetatmParser, mdk::pdb::LinkParser, mdk::pdb::ModelParser, mdk::pdb::RemarkParser, mdk::pdb::SSBondParser, mdk::pdb::TerParser

Public Functions

Record tryParse(std::string const &s) const¶

Try to parse a string into a PDB record. If unsuccessful, a std::monostate() is returned.

Parameters: s – String to parse into a PDB record.
Returns: Either a parsed PDB record or a std::monostate() if the parsing was unsuccessful.

std::ostream &write(std::ostream &os, Record const &other)¶

Write out the record (in the text form) to an output stream.

Parameters

os – Output stream to write the PDB record to.
other – PDB record to write to the stream.

Returns

os

Public Members

std::vector<std::shared_ptr<Field>> fields¶: A list of fields comprising the record.

Protected Attributes

Record record¶: Underlying record, which serves as a base of variable (integers, strings etc. of a record) to which field parsers write parsed values and from which the field parsers read the values to write.

class RemarkParser : public mdk::pdb::RecordParser ¶

Public Functions

RemarkParser()¶

class AtomParser : public mdk::pdb::RecordParser ¶

Public Functions

AtomParser()¶

class HetatmParser : public mdk::pdb::RecordParser ¶

Public Functions

HetatmParser()¶

class SSBondParser : public mdk::pdb::RecordParser ¶

Public Functions

SSBondParser()¶

class Cryst1Parser : public mdk::pdb::RecordParser ¶

Public Functions

Cryst1Parser()¶

class LinkParser : public mdk::pdb::RecordParser ¶

Public Functions

LinkParser()¶

class ModelParser : public mdk::pdb::RecordParser ¶

Public Functions

ModelParser()¶

class EndmdlParser : public mdk::pdb::RecordParser ¶

Public Functions

EndmdlParser()¶

class EndParser : public mdk::pdb::RecordParser ¶

Public Functions

EndParser()¶

class TerParser : public mdk::pdb::RecordParser ¶

Public Functions

TerParser()¶

namespace records¶

A namespace containing PODs of PDB records. For reference, see PDB format report.

Typedefs

using Record = std::variant<std::monostate, Remark, Atom, SSBond, Cryst1, End, Link, Model, Endmdl, Ter>¶: A variant type denoting either a record of some type, or no record std::monostate in the case of a parsing failure.

class Remark¶

Public Members

int number = 6¶

std::string text = ""¶

class Atom¶

Public Members

int serialNum = 0¶

std::string atomName¶

char altLocation = ' '¶

std::string residueName¶

char chainID = 'A'¶

int residueSeqNum = 0¶

char insertionCode = ' '¶

Eigen::Vector3d pos¶

double occupancy = 1.0¶

double tempFactor = 0.0¶

std::string element¶

std::string charge¶

class SSBond¶

Public Members

int serialNum = 0¶

PerResidueData res[2]¶

double dist0 = 0.0¶

struct PerResidueData¶

Public Members

char chainId = 'A'¶

int residueSeqNum = 0¶

char insertionCode = ' '¶

std::string symmetryOp¶

class Cryst1¶

Public Members

Eigen::Vector3d size¶

Eigen::Vector3d angles = {90 * degree, 90 * degree, 90 * degree}¶

std::string spaceGroup¶

int z = 0¶

class Ter¶

Public Members

int serialNum = 0¶

std::string residueName¶

char chainId = 'A'¶

int residueSeqNum = 0¶

char insertionCode = ' '¶

class Link¶

Public Members

PerResidueData res[2]¶

double linkLength¶

struct PerResidueData¶

Public Members

std::string atomName¶

char altLocation = ' '¶

std::string residueName¶

char chainId = 'A'¶

int residueSeqNum¶

char insertionCode = ' '¶

std::string symmetryOp¶

class Model¶

Public Members

int serialNum = 0¶

class Endmdl¶

class End¶

namespace seq¶

class LegacyParser¶

#include <LegacyParser.hpp>

A parser of sequence files. Here we need an explicit path since the sequence files reference contact maps by relative path.

Public Functions

Sequence read(std::filesystem::path const &path)¶

Read a sequence from a file at path.

Parameters: path – Path to a file with the sequence.
Returns: Parsed sequence.

class Sequence¶

#include <Sequence.hpp>

A sequence file, i.e. a set of chains with optional interposed contact maps.

Public Functions

Model asModel() const¶

Create a Model from the sequence file. Positions are initialized to zero, thus invoking a function like morphIntoSAW is necessary. Contacts are assumed to be “NAT”, i.e. without specifying which is sidechain or backbone etc.

Returns: Resulting Model from the sequence file

Public Members

std::optional<double> screeningDist¶: Optional screening distance, as per README.txt.

std::vector<Chain> chains¶: Chains in the sequence file.

std::unordered_map<std::string, cmap::ContactMap> contactMaps¶: A map from the contact map path to a parsed contact map.

struct Chain¶

#include <Sequence.hpp>

A single chain in the sequence file.

Public Members

std::vector<AminoAcid> codes¶: Amino acids comprising the chain.

std::vector<std::string> contactMaps¶: Paths to contact maps associated with the chain. Technically one could store cmap::ContactMap explicitly but we store paths in order to not duplicate contact maps.

namespace vl¶

struct PairInfo¶

#include <List.hpp>

A struct containing data about pairs; this is stored in this fashion chiefly to avoid having to pass all of these as arguments.

Public Members

int i1¶: Index of the first residue.

int i2¶: Index of the second residue.

Vector unit¶: Normalized vector connecting the two residues. May be non-trivial, i.e. not just r_{i_1} - r_{i_2} but the PBC-corrected version thereof.

double norm¶: Norm of the vector connecting the two residues.

class List : public mdk::SimulVar ¶

#include <List.hpp>

A Verlet list.

Public Functions

void registerNF(NonlocalForce &force, Spec const &spec)¶

Register a nonlocal force, in particular adjust the current specs so as to accomodate the newly added force (for example increase the cutoff distance), and add it to the internal list of the forces whose hooks to call when a list is reconstructed.

Parameters

force – Nonlocal force to register.
spec – Spec with which to register the force.

virtual void bind(Simulation &simulation) override¶

Bind an object to a simulation.

Parameters: simulation – Simulation to bind the object to.

void check()¶: Check whether the list needs updating, and if it does update it and invoke the relevant update hooks for the non-local forces.

Public Members

Pairs pairs¶: The list of pairs. The pairs are ordered (if a pair [i, j] is in this list, then i < j).

Private Functions

bool needToReset() const¶

void update()¶: Updates the Verlet list in the “legacy fashion”, i.e. going through a list of pairs and checking the distances.

int indexOf(Eigen::Vector3i const &loc)¶

Computes the flattened index of a cell in the grid.

Parameters: loc – Unflattened index.
Returns: Flattened index.

void perPair(int c1, int c2)¶

void perCell(int c1)¶

void updateGrid()¶: Updates the Verlet list, but instead of going through each pair we first place the residues into cells of (rough) size effCutoff and then only take the pairs from the neighboring cells, which reduces the computational cost from O(N^2) to O(N) (although the constant behind O(N) may be significant). In practice this version is much faster, especially for large numbers of residues.

Private Members

State const *state = nullptr¶

Chains const *chains = nullptr¶

double t0 = 0.0¶: Time from when the list was last reconstructed.

Vectors r0¶: Positions of the residues from when the list was last reconstructed.

Topology top0¶: Box shape from when the list was last reconstructed.

std::vector<NonlocalForce*> forces¶: A list of nonlocal forces, saved in order to invoke their update hooks when the list is reconstructed.

bool initial = false¶

double cutoff = 0.0 * angstrom ¶: A maximal cutoff among registered nonlocal forces.

double pad = 10.0 * angstrom ¶: An extra “buffer”. The list by default contains elements that may be outside the maximal cutoff, which allows us to not have to reconstruct the verlet list at every time step, at the cost of spurious distance comparisons.

int minBondSep = 0¶: Minimal bond-distance between the residues among registered nonlocal forces.

Eigen::Vector3i grid¶: Dimensions of the grid of cells used in the cell version of the computation.

std::vector<int> first¶: A list of size equal to the number of cells in the grid, where a nonnegative number indicates an index of the residue which is first in the linked list of residues in a given cell. If the cell is empty, the value is -1.

std::vector<int> last¶: A list of size equal to the number of cells in the grid, where a nonnegative number indicates an index of the residue which is last in the linked list of residues in a given cell. If the cell is empty, the value is -1.

std::vector<int> next¶: A list of size equal to the number of residues, where a nonnegative number indicates an index of a residue which is next in the linked list associated with the cell in which the residue is placed.

double effCutoff¶

double effCutoffSq¶

struct Spec¶

#include <Spec.hpp>

A specification (spec for short) for the non-local force. Technically one could generate a Verlet list for every non-local force, however this would be wasteful - thus we generate a most restrictive possible Verlet list which contains the “sub-Verlet lists” of the particular non-local forces and let them filter the general list on their own accord into a local list.

Public Members

double cutoffSq = 0.0¶: A square of the cutoff distance.

int minBondSep = 3¶: Minimum bond separation between the residues of one pair in the Verlet list. Technically this is just an optimization.