chematic_3d/pdb.rs
1//! PDB file format parser and writer.
2//!
3//! Supports `ATOM` and `HETATM` record types only. Other record types are
4//! silently ignored. Atoms are read using the standard PDB fixed-column layout
5//! (columns are 1-indexed in the PDB specification; below we use 0-indexed
6//! Rust slice indices).
7//!
8//! Column mapping (0-indexed):
9//! [0..6] record type ("ATOM " or "HETATM")
10//! [6..11] serial (right-justified integer)
11//! [12..16] atom name
12//! [17..20] residue name
13//! [21] chain ID
14//! [22..26] residue sequence number
15//! [30..38] x (8.3f)
16//! [38..46] y (8.3f)
17//! [46..54] z (8.3f)
18//! [54..60] occupancy (6.2f)
19//! [60..66] temp factor (6.2f)
20//! [76..78] element symbol
21
22use chematic_core::{Atom, AtomIdx, BondOrder, Element, Molecule, MoleculeBuilder};
23
24use crate::coords::{Coords3D, Point3};
25
26// ---------------------------------------------------------------------------
27// Public types
28// ---------------------------------------------------------------------------
29
30/// Data from a single ATOM or HETATM record in a PDB file.
31#[derive(Debug, Clone)]
32pub struct PdbAtom {
33 /// Atom serial number.
34 pub serial: u32,
35 /// Atom name field (e.g., "CA", " N ", "OXT").
36 pub name: String,
37 /// Residue name (e.g., "ALA").
38 pub res_name: String,
39 /// Chain identifier.
40 pub chain_id: char,
41 /// Residue sequence number.
42 pub res_seq: i32,
43 /// Cartesian X coordinate (angstroms).
44 pub x: f64,
45 /// Cartesian Y coordinate (angstroms).
46 pub y: f64,
47 /// Cartesian Z coordinate (angstroms).
48 pub z: f64,
49 /// Occupancy factor.
50 pub occupancy: f64,
51 /// Temperature factor (B-factor).
52 pub temp_factor: f64,
53 /// Element symbol (right-justified, 2 chars in PDB; stored trimmed here).
54 pub element: String,
55}
56
57// ---------------------------------------------------------------------------
58// Parser
59// ---------------------------------------------------------------------------
60
61/// Parse all `ATOM` and `HETATM` records from a PDB-format string.
62///
63/// Lines that are not `ATOM` / `HETATM` records (or whose record type field
64/// is not recognised) are silently skipped.
65///
66/// Atoms are capped at PDB_MAX_ATOMS (10,000) to prevent O(n²) bond inference DoS.
67pub fn parse_pdb_atoms(input: &str) -> Vec<PdbAtom> {
68 const PDB_MAX_ATOMS: usize = 10_000;
69 let mut atoms = Vec::new();
70
71 for line in input.lines() {
72 if atoms.len() >= PDB_MAX_ATOMS {
73 break;
74 }
75
76 let record = line.get(0..6).unwrap_or("").trim_end();
77 if record != "ATOM" && record != "HETATM" {
78 continue;
79 }
80
81 // Extract a fixed-width substring, tolerating short lines.
82 let field = |start: usize, end: usize| -> &str {
83 let end = end.min(line.len());
84 line.get(start..end).unwrap_or("")
85 };
86
87 atoms.push(PdbAtom {
88 serial: field(6, 11).trim().parse::<u32>().unwrap_or(0),
89 name: field(12, 16).to_string(),
90 res_name: field(17, 20).trim().to_string(),
91 chain_id: field(21, 22).chars().next().unwrap_or(' '),
92 res_seq: field(22, 26).trim().parse::<i32>().unwrap_or(0),
93 x: field(30, 38).trim().parse::<f64>().unwrap_or(0.0),
94 y: field(38, 46).trim().parse::<f64>().unwrap_or(0.0),
95 z: field(46, 54).trim().parse::<f64>().unwrap_or(0.0),
96 occupancy: field(54, 60).trim().parse::<f64>().unwrap_or(1.0),
97 temp_factor: field(60, 66).trim().parse::<f64>().unwrap_or(0.0),
98 element: field(76, 78).trim().to_string(),
99 });
100 }
101
102 atoms
103}
104
105// ---------------------------------------------------------------------------
106// PDB to Molecule
107// ---------------------------------------------------------------------------
108
109/// Build a [`Molecule`] and [`Coords3D`] from a slice of [`PdbAtom`] records.
110///
111/// Connectivity is inferred by distance: two atoms are bonded when their
112/// distance is less than `1.3 × (r_a + r_b)` where `r_a` and `r_b` are
113/// their respective covalent radii.
114pub fn pdb_to_molecule(atoms: &[PdbAtom]) -> (Molecule, Coords3D) {
115 let mut builder = MoleculeBuilder::new();
116 let mut points: Vec<Point3> = Vec::with_capacity(atoms.len());
117 let mut elements: Vec<Element> = Vec::with_capacity(atoms.len());
118
119 for pdb_atom in atoms {
120 // Use the element field when populated, else fall back to the first
121 // character of the atom name.
122 let sym = if !pdb_atom.element.is_empty() {
123 pdb_atom.element.clone()
124 } else {
125 pdb_atom
126 .name
127 .trim()
128 .chars()
129 .next()
130 .unwrap_or('C')
131 .to_string()
132 };
133 let element = Element::from_symbol(&sym).unwrap_or(Element::C);
134 elements.push(element);
135 builder.add_atom(Atom::new(element));
136 points.push(Point3::new(pdb_atom.x, pdb_atom.y, pdb_atom.z));
137 }
138
139 // Distance-based bond inference: bonded when d < 1.3 × (r_a + r_b).
140 let n = points.len();
141 for i in 0..n {
142 for j in (i + 1)..n {
143 let threshold =
144 1.3 * (elements[i].covalent_radius() as f64 + elements[j].covalent_radius() as f64);
145 if points[i].distance(&points[j]) < threshold {
146 // Ignore duplicate-bond errors (shouldn't occur for distinct pairs).
147 let _ = builder.add_bond(AtomIdx(i as u32), AtomIdx(j as u32), BondOrder::Single);
148 }
149 }
150 }
151
152 (builder.build(), Coords3D { points })
153}
154
155// ---------------------------------------------------------------------------
156// Writer
157// ---------------------------------------------------------------------------
158
159/// Write a molecule and its 3D coordinates as a PDB string (HETATM records).
160///
161/// All atoms are written as `HETATM` records in a single residue ("LIG", chain A,
162/// residue 1). No `CONECT` records are written.
163///
164/// Column layout (0-indexed):
165/// [0..6] "HETATM"
166/// [6..11] serial (5 chars, right-justified)
167/// [11] space
168/// [12..16] atom name (4 chars, e.g. " C ")
169/// [16] altLoc (space)
170/// [17..20] resName (3 chars)
171/// [20] space
172/// [21] chainID
173/// [22..26] resSeq (4 chars, right-justified)
174/// [26] iCode (space)
175/// [27..30] 3 spaces
176/// [30..38] x (8.3f)
177/// [38..46] y (8.3f)
178/// [46..54] z (8.3f)
179/// [54..60] occupancy (6.2f)
180/// [60..66] tempFactor (6.2f)
181/// [66..76] 10 spaces
182/// [76..78] element (2 chars, right-justified)
183pub fn write_pdb(mol: &Molecule, coords: &Coords3D) -> String {
184 let mut out = String::new();
185
186 for i in 0..mol.atom_count() {
187 let idx = AtomIdx(i as u32);
188 let atom = mol.atom(idx);
189 let p = coords.get(idx);
190
191 let serial = (i + 1) as u32;
192 let elem_sym = atom.element.symbol();
193
194 // PDB atom name field (4 chars): single-char elements left-padded with space.
195 let atom_name = if elem_sym.len() == 1 {
196 format!(" {:<3}", elem_sym)
197 } else {
198 format!("{:<4}", elem_sym)
199 };
200
201 // Build the line by assembling fixed-width fields directly.
202 // Total width must reach at least col 78 for the element field.
203 let line = format!(
204 // 0-5 6-10 11 12-15 16 17-19 20 21 22-25 26 27-29 30-37 38-45 46-53 54-59 60-65 66-75 76-77
205 "{:<6}{:>5} {:<4} {:>3} {:1}{:>4} {:>8.3}{:>8.3}{:>8.3}{:>6.2}{:>6.2} {:>2}\n",
206 "HETATM", // 0-5
207 serial, // 6-10
208 // space at 11 (literal in format)
209 atom_name, // 12-15
210 // altLoc at 16 (space, literal in format)
211 "LIG", // 17-19
212 // space at 20 (literal in format)
213 'A', // 21
214 1_i32, // 22-25
215 // iCode + 3 spaces (27-29) handled by 4-space literal in format
216 p.x, // 30-37
217 p.y, // 38-45
218 p.z, // 46-53
219 1.00_f64, // 54-59
220 0.00_f64, // 60-65
221 elem_sym, // 76-77
222 );
223 out.push_str(&line);
224 }
225
226 out.push_str("END\n");
227 out
228}