Generate elastic network ======================== The ``genelastic`` tool generates **elastic network bonds** for coarse-grained simulation topologies (``.itp`` files). Elastic networks are widely used in CG/HPS/Go-like models to enforce: - tertiary structure stability - domain rigidity - elastic network modeling (ENM) - local structural preservation across residues ``genelastic`` reads: - a **PDB** structure (for atom coordinates), - an input **ITP** topology (for atom list & existing bonds), - an **elastic residue cluster file** (ASCII list of atom groups), and generates new *Bond* entries between atom pairs whose distances fall within user-specified bounds. This tool is implemented in ``genelastic.py``. Overview -------- The program builds an elastic network by: 1. Loading atomic coordinates from a **PDB** file 2. Loading a topology from an **ITP** file 3. Reading an ASCII file specifying **clusters of residue IDs** 4. For each cluster: - consider all atom pairs - skip pairs already bonded in topology - skip pairs already added earlier in this run - compute the distance from PDB coordinates - if the distance is within the cutoff window ``elastic_lower ≤ r ≤ elastic_upper`` → add an elastic bond 5. Write updated topology to an output ITP file Elastic bonds are harmonic: .. math:: V(r) = \frac{1}{2} k(r - r_0)^2 where ``k`` is the elastic-force constant and ``r_0`` is the native contact distance. Usage ----- Minimal example: .. code-block:: bash dps genelastic \ -f structure.pdb \ -p input.itp \ -o output.itp \ -er clusters.txt Arguments --------- Required arguments ~~~~~~~~~~~~~~~~~~ .. option:: -f, --structure FILE Input PDB structure. Coordinates are used to determine native distances. .. option:: -p, --topology FILE Input ITP file whose topology is to be extended. Must end with ``.itp``. .. option:: -o, --output FILE Output ITP filename. Extension ``.itp`` is added automatically if missing. .. option:: -er, --elastic-residues FILE ASCII file where each line lists a **cluster** of atoms (1-based indices). Example: :: 3 5 8 10 2 4 7 Elastic-network parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~ .. option:: -ef, --elastic-force-constant FLOAT Force constant (``k``) in **kJ / mol / nm²**. Default: ``5000``. .. option:: -el, --elastic-lower FLOAT Lower distance cutoff in **nm**. Only pairs with ``r > elastic_lower`` are considered. Default: ``0.5``. .. option:: -eu, --elastic-upper FLOAT Upper distance cutoff in **nm**. Only pairs with ``r < elastic_upper`` are considered. Default: ``0.9``. Input File Formats ------------------ ### 1. PDB File (structure) Used only for coordinates (x, y, z) and box vectors. ### 2. ITP File (topology) Must contain: - ``[ atoms ]`` - ``[ bonds ]`` The elastic network is added to the ``[ bonds ]`` section using: .. code-block:: python Bond(atom1, atom2, bond_length, force_constant) ### 3. Elastic Cluster File (ASCII) Each non-empty line: :: idx1 idx2 idx3 ... These indexes are 1-based **atom indices**, *not residue numbers*. A **cluster** means *every atom in the line can potentially form an elastic bond with every other atom on that line*, subject to distance criteria. Program Workflow ---------------- (1) **Read PDB structure** .. code-block:: python atoms, box = read_pdb(args.structure) If reading fails: :: ## An exception occurred when trying to open structure file... (2) **Read ITP topology** .. code-block:: python topology = read_itp(args.topology) If input ITP file does not end in ``.itp`` → immediate error. (3) **Read elastic clusters** .. code-block:: python elastic_cluster_list = [ [int(i)-1 for i in line.split()] for line in open(args.elastic_residues) ] Each atom index is converted to **0-based**. (4) **Prepare existing bond lookup** .. code-block:: python exist_bonds = [[bond.a1, bond.a2] for bond in topology.bonds] Used to avoid overwriting or duplicating bonds. (5) **Loop over each cluster & atom pair** For every pair ``(i, j)`` in a cluster: - skip if ``i == j`` - skip if pair appears in ``exist_bonds`` - skip if pair appears in previously ``added_bonds`` - compute distance from PDB coordinates: .. code-block:: python distance = np.linalg.norm(coor_1 - coor_2) - check: :: elastic_lower ≤ distance ≤ elastic_upper - if passed → add an elastic bond with: - ``bond_length = distance`` (in nm) - ``k`` = elastic-force-constant (converted to OpenMM units) Then: .. code-block:: python topology.bonds.append(Bond(atomid_1, atomid_2, bond_length, bond_k)) Program message: :: ## Adding bond (k=...) between atom X and Y with distance Z. (6) **Save output** .. code-block:: python write_itp(output_file_name, topology) Final message: :: ## A total number of N bonds have been added. ## Output topology file written to .itp. Example ------- Example elastic residue file: :: 5 6 7 10 20 30 35 Run: .. code-block:: bash dps genelastic \ -f model.pdb \ -p protein.itp \ -o protein_elastic.itp \ -er elastic_groups.txt \ -el 0.4 -eu 1.0 -ef 6000 Example output: :: ## Open structure file model.pdb which contains 414 atoms. ## Open topology file protein.itp which contains 414 atoms. ## Reference structure contains 2 clusters. ## Will generate elastic bonds between residue within 0.40 ~ 1.00 nm. ## Elastic bond constant will be 6000.00 kJ / nm ^ 2 ## Adding bond (k=...) between atom 5 and 6 with distance 0.87. ## Adding bond (k=...) between atom 5 and 7 with distance 0.71. ... ## A total number of 48 bonds have been added. ## Output topology file written to protein_elastic.itp. Error Messages -------------- **“ERROR: Only itp file can be treated as input.”** Input topology must be ``.itp`` format. **“An exception occurred when trying to open structure file …”** PDB could not be loaded. **“An exception occurred when trying to open angle file …”** Actually refers to elastic-residue file loading error. **Duplicate or conflicting bonds** Automatically skipped silently (existing bonds take priority). Summary ------- ``dps genelastic`` is a powerful generator of elastic networks for coarse-grained simulations. It provides: - cluster-based ENM construction - distance-based filtering - duplicate-bond protection - clean integration into existing ITP files - precise unit handling (nm, kJ/mol/nm²) It is ideal for stabilizing tertiary structures or imposing distance constraints in CGPS / DROPPS simulations.