zsasa
zsasa: Fast SASA calculation using Zig.
This package provides Python bindings for the zsasa library, a high-performance implementation of Solvent Accessible Surface Area (SASA) calculation algorithms.
Example:
import numpy as np from zsasa import calculate_sasa
Single atom
coords = np.array([[0.0, 0.0, 0.0]]) radii = np.array([1.5]) result = calculate_sasa(coords, radii) print(f"Total SASA: {result.total_area:.2f} Ų")
>>> # Classify atoms >>> from zsasa import classify_atoms, get_radius >>> result = classify_atoms(["ALA", "ALA"], ["CA", "O"]) >>> print(result.radii) # [1.88 1.42] (default CCD classifier)Integrations: For structure file support, use the gemmi integration:
>>> # pip install zsasa[gemmi] >>> from zsasa.integrations.gemmi import calculate_sasa_from_structure >>> result = calculate_sasa_from_structure("protein.cif") >>> print(f"Total: {result.total_area:.1f} Ų")Analysis: For per-residue aggregation and RSA calculation:
>>> from zsasa import aggregate_from_result >>> from zsasa.integrations.gemmi import calculate_sasa_from_structure >>> result = calculate_sasa_from_structure("protein.cif") >>> residues = aggregate_from_result(result) >>> for res in residues: ... if res.rsa is not None: ... print(f"{res.chain_id}:{res.residue_name}{res.residue_id}: RSA={res.rsa:.1%}")MDTraj Integration: For MD trajectory analysis (requires mdtraj):
>>> # pip install mdtraj >>> from zsasa.mdtraj import compute_sasa >>> import mdtraj as md >>> traj = md.load('trajectory.xtc', top='topology.pdb') >>> sasa = compute_sasa(traj) # Returns (n_frames, n_atoms) in nm²MDAnalysis Integration: For MD trajectory analysis with MDAnalysis (requires MDAnalysis):
>>> # pip install MDAnalysis >>> import MDAnalysis as mda >>> from zsasa.mdanalysis import SASAAnalysis >>> u = mda.Universe('topology.pdb', 'trajectory.xtc') >>> sasa = SASAAnalysis(u, select='protein') >>> sasa.run() >>> print(sasa.results.total_area) # Returns per-frame SASA in ŲTrajectory File I/O: For direct trajectory-file I/O and trajectory-native analysis in Python, prefer pyztraj (the Python bindings for ztraj). zsasa keeps the legacy zsasa.xtc and zsasa.dcd modules for compatibility, but new trajectory formats are centralized in pyztraj.
>>> # pip install pyztraj >>> import pyztraj >>> structure = pyztraj.load_pdb("topology.pdb") >>> with pyztraj.open_trr("trajectory.trr", structure.n_atoms) as reader: ... for frame in reader: ... sasa = pyztraj.compute_sasa(structure, frame.coords) ... print(sasa.total_area)
1"""zsasa: Fast SASA calculation using Zig. 2 3This package provides Python bindings for the zsasa library, 4a high-performance implementation of Solvent Accessible Surface Area (SASA) 5calculation algorithms. 6 7Example: 8 >>> import numpy as np 9 >>> from zsasa import calculate_sasa 10 >>> 11 >>> # Single atom 12 >>> coords = np.array([[0.0, 0.0, 0.0]]) 13 >>> radii = np.array([1.5]) 14 >>> result = calculate_sasa(coords, radii) 15 >>> print(f"Total SASA: {result.total_area:.2f} Ų") 16 17 >>> # Classify atoms 18 >>> from zsasa import classify_atoms, get_radius 19 >>> result = classify_atoms(["ALA", "ALA"], ["CA", "O"]) 20 >>> print(result.radii) # [1.88 1.42] (default CCD classifier) 21 22Integrations: 23 For structure file support, use the gemmi integration: 24 25 >>> # pip install zsasa[gemmi] 26 >>> from zsasa.integrations.gemmi import calculate_sasa_from_structure 27 >>> result = calculate_sasa_from_structure("protein.cif") 28 >>> print(f"Total: {result.total_area:.1f} Ų") 29 30Analysis: 31 For per-residue aggregation and RSA calculation: 32 33 >>> from zsasa import aggregate_from_result 34 >>> from zsasa.integrations.gemmi import calculate_sasa_from_structure 35 >>> result = calculate_sasa_from_structure("protein.cif") 36 >>> residues = aggregate_from_result(result) 37 >>> for res in residues: 38 ... if res.rsa is not None: 39 ... print(f"{res.chain_id}:{res.residue_name}{res.residue_id}: RSA={res.rsa:.1%}") 40 41MDTraj Integration: 42 For MD trajectory analysis (requires mdtraj): 43 44 >>> # pip install mdtraj 45 >>> from zsasa.mdtraj import compute_sasa 46 >>> import mdtraj as md 47 >>> traj = md.load('trajectory.xtc', top='topology.pdb') 48 >>> sasa = compute_sasa(traj) # Returns (n_frames, n_atoms) in nm² 49 50MDAnalysis Integration: 51 For MD trajectory analysis with MDAnalysis (requires MDAnalysis): 52 53 >>> # pip install MDAnalysis 54 >>> import MDAnalysis as mda 55 >>> from zsasa.mdanalysis import SASAAnalysis 56 >>> u = mda.Universe('topology.pdb', 'trajectory.xtc') 57 >>> sasa = SASAAnalysis(u, select='protein') 58 >>> sasa.run() 59 >>> print(sasa.results.total_area) # Returns per-frame SASA in Ų 60 61Trajectory File I/O: 62 For direct trajectory-file I/O and trajectory-native analysis in Python, 63 prefer pyztraj (the Python bindings for ztraj). zsasa keeps the legacy 64 zsasa.xtc and zsasa.dcd modules for compatibility, but new trajectory 65 formats are centralized in pyztraj. 66 67 >>> # pip install pyztraj 68 >>> import pyztraj 69 >>> structure = pyztraj.load_pdb("topology.pdb") 70 >>> with pyztraj.open_trr("trajectory.trr", structure.n_atoms) as reader: 71 ... for frame in reader: 72 ... sasa = pyztraj.compute_sasa(structure, frame.coords) 73 ... print(sasa.total_area) 74""" 75 76from zsasa._ffi import get_version 77from zsasa.analysis import ( 78 ResidueResult, 79 aggregate_by_residue, 80 aggregate_from_result, 81) 82from zsasa.batch import BatchDirResult, process_directory 83from zsasa.classifier import ( 84 AtomClass, 85 ClassificationResult, 86 ClassifierType, 87 classify_atoms, 88 get_atom_class, 89 get_radius, 90 guess_radius, 91 guess_radius_from_atom_name, 92) 93from zsasa.rsa import ( 94 MAX_SASA, 95 calculate_rsa, 96 calculate_rsa_batch, 97 get_max_sasa, 98) 99from zsasa.sasa import ( 100 BatchSasaResult, 101 SasaResult, 102 calculate_sasa, 103 calculate_sasa_batch, 104) 105 106__all__ = [ 107 # SASA calculation 108 "calculate_sasa", 109 "calculate_sasa_batch", 110 "SasaResult", 111 "BatchSasaResult", 112 # Batch directory processing 113 "process_directory", 114 "BatchDirResult", 115 # Classifier 116 "ClassifierType", 117 "AtomClass", 118 "ClassificationResult", 119 "get_radius", 120 "get_atom_class", 121 "guess_radius", 122 "guess_radius_from_atom_name", 123 "classify_atoms", 124 # RSA 125 "MAX_SASA", 126 "get_max_sasa", 127 "calculate_rsa", 128 "calculate_rsa_batch", 129 # Analysis 130 "ResidueResult", 131 "aggregate_by_residue", 132 "aggregate_from_result", 133 # Utility 134 "get_version", 135] 136 137__version__ = get_version()
106def calculate_sasa( 107 coords: NDArray[np.float64], 108 radii: NDArray[np.float64], 109 *, 110 algorithm: Literal["sr", "lr"] = "sr", 111 n_points: int = 100, 112 n_slices: int = 20, 113 probe_radius: float = 1.4, 114 n_threads: int = 0, 115 use_bitmask: bool = False, 116 bitmask_correction: bool = False, 117 bitmask_correction_coeff: float | None = None, 118) -> SasaResult: 119 """Calculate Solvent Accessible Surface Area (SASA). 120 121 Args: 122 coords: Atom coordinates as (N, 3) array. 123 radii: Atom radii as (N,) array. 124 algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards). 125 Lee-Richards computes its arc angles with exact trigonometry; the 126 approximate mode of the command line (--lr-trig=fast) is not 127 available from Python. 128 n_points: Number of test points per atom (for SR algorithm). Default: 100. 129 n_slices: Number of slices per atom (for LR algorithm). Default: 20. 130 probe_radius: Water probe radius in Angstroms. Default: 1.4. 131 n_threads: Number of threads to use. 0 = auto-detect. Default: 0. 132 use_bitmask: Use bitmask LUT optimization for SR algorithm. 133 Supports n_points 1..1024. Default: False. 134 bitmask_correction: Apply experimental bitmask exposed-fraction 135 correction. Requires use_bitmask=True. Default: False. 136 bitmask_correction_coeff: Optional non-negative correction coefficient. 137 Defaults to the library's experimental coefficient when omitted. 138 139 Returns: 140 SasaResult containing total_area and per-atom atom_areas. 141 142 Raises: 143 ValueError: If input arrays have invalid shapes or calculation fails. 144 145 Example: 146 >>> import numpy as np 147 >>> from zsasa import calculate_sasa 148 >>> 149 >>> # Two atoms 150 >>> coords = np.array([[0.0, 0.0, 0.0], [3.0, 0.0, 0.0]]) 151 >>> radii = np.array([1.5, 1.5]) 152 >>> result = calculate_sasa(coords, radii) 153 >>> print(f"Total: {result.total_area:.2f}") 154 """ 155 # Validate parameters before loading the C library so invalid Python inputs 156 # never reach CFFI conversion or native code. 157 n_points = _validate_uint32_param("n_points", n_points) 158 n_slices = _validate_uint32_param("n_slices", n_slices) 159 probe_radius = _validate_positive_finite_float("probe_radius", probe_radius) 160 n_threads = _validate_size_t_param("n_threads", n_threads) 161 if bitmask_correction and not use_bitmask: 162 msg = "bitmask_correction=True requires use_bitmask=True" 163 raise ValueError(msg) 164 if bitmask_correction_coeff is not None and ( 165 bitmask_correction_coeff < 0 or not np.isfinite(bitmask_correction_coeff) 166 ): 167 msg = f"bitmask_correction_coeff must be non-negative, got {bitmask_correction_coeff}" 168 raise ValueError(msg) 169 170 # Validate bitmask constraints 171 if use_bitmask: 172 use_bitmask = _validate_bitmask_params(algorithm, n_points, strict=bitmask_correction) 173 174 # Validate and convert inputs 175 coords = np.ascontiguousarray(coords, dtype=np.float64) 176 radii = np.ascontiguousarray(radii, dtype=np.float64) 177 178 if coords.ndim != 2 or coords.shape[1] != 3: 179 msg = f"coords must be (N, 3) array, got shape {coords.shape}" 180 raise ValueError(msg) 181 182 n_atoms = coords.shape[0] 183 if radii.shape != (n_atoms,): 184 msg = f"radii must be ({n_atoms},) array, got shape {radii.shape}" 185 raise ValueError(msg) 186 187 _validate_finite_array("coords", coords) 188 _validate_finite_array("radii", radii) 189 190 if np.any(radii < 0): 191 msg = "All radii must be non-negative" 192 raise ValueError(msg) 193 194 ffi, lib = _get_lib() 195 196 # Extract x, y, z as contiguous arrays 197 x = np.ascontiguousarray(coords[:, 0]) 198 y = np.ascontiguousarray(coords[:, 1]) 199 z = np.ascontiguousarray(coords[:, 2]) 200 201 # Allocate output arrays 202 atom_areas = np.zeros(n_atoms, dtype=np.float64) 203 total_area = ffi.new("double*") 204 205 # Get cffi pointers from numpy arrays 206 x_ptr = ffi.cast("double*", x.ctypes.data) 207 y_ptr = ffi.cast("double*", y.ctypes.data) 208 z_ptr = ffi.cast("double*", z.ctypes.data) 209 radii_ptr = ffi.cast("double*", radii.ctypes.data) 210 areas_ptr = ffi.cast("double*", atom_areas.ctypes.data) 211 212 # Call the appropriate function 213 if use_bitmask: 214 if bitmask_correction: 215 coeff = 0.020 if bitmask_correction_coeff is None else bitmask_correction_coeff 216 result = lib.zsasa_calc_sr_bitmask_corrected( 217 x_ptr, 218 y_ptr, 219 z_ptr, 220 radii_ptr, 221 n_atoms, 222 n_points, 223 probe_radius, 224 n_threads, 225 coeff, 226 areas_ptr, 227 total_area, 228 ) 229 else: 230 result = lib.zsasa_calc_sr_bitmask( 231 x_ptr, 232 y_ptr, 233 z_ptr, 234 radii_ptr, 235 n_atoms, 236 n_points, 237 probe_radius, 238 n_threads, 239 areas_ptr, 240 total_area, 241 ) 242 elif algorithm == "sr": 243 result = lib.zsasa_calc_sr( 244 x_ptr, 245 y_ptr, 246 z_ptr, 247 radii_ptr, 248 n_atoms, 249 n_points, 250 probe_radius, 251 n_threads, 252 areas_ptr, 253 total_area, 254 ) 255 elif algorithm == "lr": 256 result = lib.zsasa_calc_lr( 257 x_ptr, 258 y_ptr, 259 z_ptr, 260 radii_ptr, 261 n_atoms, 262 n_slices, 263 probe_radius, 264 n_threads, 265 areas_ptr, 266 total_area, 267 ) 268 else: 269 msg = f"Unknown algorithm: {algorithm}. Use 'sr' or 'lr'." 270 raise ValueError(msg) 271 272 # Check for errors 273 if result == ZSASA_ERROR_INVALID_INPUT: 274 msg = "Invalid input to SASA calculation" 275 raise ValueError(msg) 276 elif result == ZSASA_ERROR_OUT_OF_MEMORY: 277 msg = "Out of memory during SASA calculation" 278 raise MemoryError(msg) 279 elif result == ZSASA_ERROR_CALCULATION: 280 msg = "Error during SASA calculation" 281 raise RuntimeError(msg) 282 elif result == ZSASA_ERROR_UNSUPPORTED_N_POINTS: 283 msg = ( 284 f"Unsupported n_points for bitmask: {n_points}. " 285 f"Must be {_BITMASK_MIN_N_POINTS}..{_BITMASK_MAX_N_POINTS}" 286 ) 287 raise ValueError(msg) 288 elif result != ZSASA_OK: 289 msg = f"Unknown error code: {result}" 290 raise RuntimeError(msg) 291 292 return SasaResult( 293 total_area=total_area[0], 294 atom_areas=atom_areas, 295 )
Calculate Solvent Accessible Surface Area (SASA).
Args: coords: Atom coordinates as (N, 3) array. radii: Atom radii as (N,) array. algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards). Lee-Richards computes its arc angles with exact trigonometry; the approximate mode of the command line (--lr-trig=fast) is not available from Python. n_points: Number of test points per atom (for SR algorithm). Default: 100. n_slices: Number of slices per atom (for LR algorithm). Default: 20. probe_radius: Water probe radius in Angstroms. Default: 1.4. n_threads: Number of threads to use. 0 = auto-detect. Default: 0. use_bitmask: Use bitmask LUT optimization for SR algorithm. Supports n_points 1..1024. Default: False. bitmask_correction: Apply experimental bitmask exposed-fraction correction. Requires use_bitmask=True. Default: False. bitmask_correction_coeff: Optional non-negative correction coefficient. Defaults to the library's experimental coefficient when omitted.
Returns: SasaResult containing total_area and per-atom atom_areas.
Raises: ValueError: If input arrays have invalid shapes or calculation fails.
Example:
import numpy as np from zsasa import calculate_sasa
Two atoms
coords = np.array([[0.0, 0.0, 0.0], [3.0, 0.0, 0.0]]) radii = np.array([1.5, 1.5]) result = calculate_sasa(coords, radii) print(f"Total: {result.total_area:.2f}")
328def calculate_sasa_batch( 329 coordinates: NDArray[np.floating], 330 radii: NDArray[np.floating], 331 *, 332 algorithm: Literal["sr", "lr"] = "sr", 333 n_points: int = 100, 334 n_slices: int = 20, 335 probe_radius: float = 1.4, 336 n_threads: int = 0, 337 precision: Literal["f64", "f32"] = "f64", 338 use_bitmask: bool = False, 339 bitmask_correction: bool = False, 340 bitmask_correction_coeff: float | None = None, 341) -> BatchSasaResult: 342 """Calculate SASA for multiple frames (batch processing). 343 344 Optimized for MD trajectory analysis where the same atoms are processed 345 across multiple frames. Parallelizes across frames for maximum performance. 346 347 Args: 348 coordinates: Atom coordinates as (n_frames, n_atoms, 3) array in Angstroms. 349 radii: Atom radii as (n_atoms,) array in Angstroms. 350 algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards). 351 n_points: Number of test points per atom (for SR algorithm). Default: 100. 352 n_slices: Number of slices per atom (for LR algorithm). Default: 20. 353 probe_radius: Water probe radius in Angstroms. Default: 1.4. 354 n_threads: Number of threads to use. 0 = auto-detect. Default: 0. 355 precision: Internal calculation precision: "f64" (default, higher precision) 356 or "f32" (matches RustSASA/mdsasa-bolt for comparison). Default: "f64". 357 use_bitmask: Use bitmask LUT optimization for SR algorithm. 358 Supports n_points 1..1024. Default: False. 359 bitmask_correction: Apply experimental bitmask exposed-fraction 360 correction. Requires use_bitmask=True. Default: False. 361 bitmask_correction_coeff: Optional non-negative correction coefficient. 362 Defaults to the library's experimental coefficient when omitted. 363 364 Returns: 365 BatchSasaResult containing per-atom SASA for all frames. 366 367 Raises: 368 ValueError: If input arrays have invalid shapes or calculation fails. 369 370 Example: 371 >>> import numpy as np 372 >>> from zsasa import calculate_sasa_batch 373 >>> 374 >>> # 10 frames, 100 atoms 375 >>> coords = np.random.randn(10, 100, 3).astype(np.float32) 376 >>> radii = np.full(100, 1.5, dtype=np.float32) 377 >>> result = calculate_sasa_batch(coords, radii) 378 >>> print(f"Shape: {result.atom_areas.shape}") # (10, 100) 379 >>> print(f"Total SASA per frame: {result.total_areas}") 380 """ 381 # Validate parameters before loading the C library so invalid Python inputs 382 # never reach CFFI conversion or native code. 383 n_points = _validate_uint32_param("n_points", n_points) 384 n_slices = _validate_uint32_param("n_slices", n_slices) 385 probe_radius = _validate_positive_finite_float("probe_radius", probe_radius) 386 n_threads = _validate_size_t_param("n_threads", n_threads) 387 if bitmask_correction and not use_bitmask: 388 msg = "bitmask_correction=True requires use_bitmask=True" 389 raise ValueError(msg) 390 if bitmask_correction_coeff is not None and ( 391 bitmask_correction_coeff < 0 or not np.isfinite(bitmask_correction_coeff) 392 ): 393 msg = f"bitmask_correction_coeff must be non-negative, got {bitmask_correction_coeff}" 394 raise ValueError(msg) 395 396 # Validate bitmask constraints 397 if use_bitmask: 398 use_bitmask = _validate_bitmask_params(algorithm, n_points, strict=bitmask_correction) 399 400 # Validate and convert inputs 401 coordinates = np.ascontiguousarray(coordinates, dtype=np.float32) 402 radii = np.ascontiguousarray(radii, dtype=np.float32) 403 404 if coordinates.ndim != 3 or coordinates.shape[2] != 3: 405 msg = f"coordinates must be (n_frames, n_atoms, 3) array, got shape {coordinates.shape}" 406 raise ValueError(msg) 407 408 n_frames = coordinates.shape[0] 409 n_atoms = coordinates.shape[1] 410 411 if radii.shape != (n_atoms,): 412 msg = f"radii must be ({n_atoms},) array, got shape {radii.shape}" 413 raise ValueError(msg) 414 415 _validate_finite_array("coordinates", coordinates) 416 _validate_finite_array("radii", radii) 417 418 if np.any(radii < 0): 419 msg = "All radii must be non-negative" 420 raise ValueError(msg) 421 422 ffi, lib = _get_lib() 423 424 # Allocate output array 425 atom_areas = np.zeros((n_frames, n_atoms), dtype=np.float32) 426 427 # Get cffi pointers from numpy arrays 428 coords_ptr = ffi.cast("float*", coordinates.ctypes.data) 429 radii_ptr = ffi.cast("float*", radii.ctypes.data) 430 areas_ptr = ffi.cast("float*", atom_areas.ctypes.data) 431 432 # Call the appropriate batch function based on algorithm, precision, and bitmask 433 if use_bitmask: 434 coeff = 0.020 if bitmask_correction_coeff is None else bitmask_correction_coeff 435 if precision == "f32": 436 if bitmask_correction: 437 result = lib.zsasa_calc_sr_batch_bitmask_f32_corrected( 438 coords_ptr, 439 n_frames, 440 n_atoms, 441 radii_ptr, 442 n_points, 443 probe_radius, 444 n_threads, 445 coeff, 446 areas_ptr, 447 ) 448 else: 449 result = lib.zsasa_calc_sr_batch_bitmask_f32( 450 coords_ptr, 451 n_frames, 452 n_atoms, 453 radii_ptr, 454 n_points, 455 probe_radius, 456 n_threads, 457 areas_ptr, 458 ) 459 else: 460 if bitmask_correction: 461 result = lib.zsasa_calc_sr_batch_bitmask_corrected( 462 coords_ptr, 463 n_frames, 464 n_atoms, 465 radii_ptr, 466 n_points, 467 probe_radius, 468 n_threads, 469 coeff, 470 areas_ptr, 471 ) 472 else: 473 result = lib.zsasa_calc_sr_batch_bitmask( 474 coords_ptr, 475 n_frames, 476 n_atoms, 477 radii_ptr, 478 n_points, 479 probe_radius, 480 n_threads, 481 areas_ptr, 482 ) 483 elif precision == "f64": 484 # Default: f32 I/O with f64 internal precision 485 if algorithm == "sr": 486 result = lib.zsasa_calc_sr_batch( 487 coords_ptr, 488 n_frames, 489 n_atoms, 490 radii_ptr, 491 n_points, 492 probe_radius, 493 n_threads, 494 areas_ptr, 495 ) 496 elif algorithm == "lr": 497 result = lib.zsasa_calc_lr_batch( 498 coords_ptr, 499 n_frames, 500 n_atoms, 501 radii_ptr, 502 n_slices, 503 probe_radius, 504 n_threads, 505 areas_ptr, 506 ) 507 else: 508 msg = f"Unknown algorithm: {algorithm}. Use 'sr' or 'lr'." 509 raise ValueError(msg) 510 elif precision == "f32": 511 # Pure f32 precision (matches RustSASA/mdsasa-bolt) 512 if algorithm == "sr": 513 result = lib.zsasa_calc_sr_batch_f32( 514 coords_ptr, 515 n_frames, 516 n_atoms, 517 radii_ptr, 518 n_points, 519 probe_radius, 520 n_threads, 521 areas_ptr, 522 ) 523 elif algorithm == "lr": 524 result = lib.zsasa_calc_lr_batch_f32( 525 coords_ptr, 526 n_frames, 527 n_atoms, 528 radii_ptr, 529 n_slices, 530 probe_radius, 531 n_threads, 532 areas_ptr, 533 ) 534 else: 535 msg = f"Unknown algorithm: {algorithm}. Use 'sr' or 'lr'." 536 raise ValueError(msg) 537 else: 538 msg = f"Unknown precision: {precision}. Use 'f64' or 'f32'." 539 raise ValueError(msg) 540 541 # Check for errors 542 if result == ZSASA_ERROR_INVALID_INPUT: 543 msg = "Invalid input to batch SASA calculation" 544 raise ValueError(msg) 545 elif result == ZSASA_ERROR_OUT_OF_MEMORY: 546 msg = "Out of memory during batch SASA calculation" 547 raise MemoryError(msg) 548 elif result == ZSASA_ERROR_CALCULATION: 549 msg = "Error during batch SASA calculation" 550 raise RuntimeError(msg) 551 elif result == ZSASA_ERROR_UNSUPPORTED_N_POINTS: 552 msg = ( 553 f"Unsupported n_points for bitmask: {n_points}. " 554 f"Must be {_BITMASK_MIN_N_POINTS}..{_BITMASK_MAX_N_POINTS}" 555 ) 556 raise ValueError(msg) 557 elif result != ZSASA_OK: 558 msg = f"Unknown error code: {result}" 559 raise RuntimeError(msg) 560 561 return BatchSasaResult(atom_areas=atom_areas)
Calculate SASA for multiple frames (batch processing).
Optimized for MD trajectory analysis where the same atoms are processed across multiple frames. Parallelizes across frames for maximum performance.
Args: coordinates: Atom coordinates as (n_frames, n_atoms, 3) array in Angstroms. radii: Atom radii as (n_atoms,) array in Angstroms. algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards). n_points: Number of test points per atom (for SR algorithm). Default: 100. n_slices: Number of slices per atom (for LR algorithm). Default: 20. probe_radius: Water probe radius in Angstroms. Default: 1.4. n_threads: Number of threads to use. 0 = auto-detect. Default: 0. precision: Internal calculation precision: "f64" (default, higher precision) or "f32" (matches RustSASA/mdsasa-bolt for comparison). Default: "f64". use_bitmask: Use bitmask LUT optimization for SR algorithm. Supports n_points 1..1024. Default: False. bitmask_correction: Apply experimental bitmask exposed-fraction correction. Requires use_bitmask=True. Default: False. bitmask_correction_coeff: Optional non-negative correction coefficient. Defaults to the library's experimental coefficient when omitted.
Returns: BatchSasaResult containing per-atom SASA for all frames.
Raises: ValueError: If input arrays have invalid shapes or calculation fails.
Example:
import numpy as np from zsasa import calculate_sasa_batch
10 frames, 100 atoms
coords = np.random.randn(10, 100, 3).astype(np.float32) radii = np.full(100, 1.5, dtype=np.float32) result = calculate_sasa_batch(coords, radii) print(f"Shape: {result.atom_areas.shape}") # (10, 100) print(f"Total SASA per frame: {result.total_areas}")
93@dataclass 94class SasaResult: 95 """Result of SASA calculation. 96 97 Attributes: 98 total_area: Total solvent accessible surface area in Ų. 99 atom_areas: Per-atom SASA values in Ų. 100 """ 101 102 total_area: float 103 atom_areas: NDArray[np.float64]
Result of SASA calculation.
Attributes: total_area: Total solvent accessible surface area in Ų. atom_areas: Per-atom SASA values in Ų.
298@dataclass 299class BatchSasaResult: 300 """Result of batch SASA calculation for multiple frames. 301 302 Attributes: 303 atom_areas: Per-atom SASA values for all frames, shape (n_frames, n_atoms). 304 Values are in Angstrom² (Ų). 305 """ 306 307 atom_areas: NDArray[np.float32] 308 309 @property 310 def n_frames(self) -> int: 311 """Number of frames.""" 312 return self.atom_areas.shape[0] 313 314 @property 315 def n_atoms(self) -> int: 316 """Number of atoms.""" 317 return self.atom_areas.shape[1] 318 319 @property 320 def total_areas(self) -> NDArray[np.float32]: 321 """Total SASA per frame, shape (n_frames,).""" 322 return self.atom_areas.sum(axis=1) 323 324 def __repr__(self) -> str: 325 return f"BatchSasaResult(n_frames={self.n_frames}, n_atoms={self.n_atoms})"
Result of batch SASA calculation for multiple frames.
Attributes: atom_areas: Per-atom SASA values for all frames, shape (n_frames, n_atoms). Values are in Angstrom² (Ų).
309 @property 310 def n_frames(self) -> int: 311 """Number of frames.""" 312 return self.atom_areas.shape[0]
Number of frames.
109def process_directory( 110 input_dir: str | Path, 111 *, 112 output_dir: str | Path | None = None, 113 algorithm: Literal["sr", "lr"] = "sr", 114 n_points: int = 100, 115 n_slices: int = 20, 116 probe_radius: float = 1.4, 117 n_threads: int = 0, 118 classifier: ClassifierType | None = ClassifierType.CCD, 119 include_hydrogens: bool = False, 120 include_hetatm: bool = False, 121) -> BatchDirResult: 122 """Process all supported structure files in a directory for SASA calculation. 123 124 Supported formats: PDB (.pdb, .ent), mmCIF (.cif, .mmcif), BinaryCIF (.bcif), 125 JSON (.json), SDF/MOL (.sdf, .mol), and their gzip- or zstd-compressed 126 variants (.gz, .zst). 127 128 Args: 129 input_dir: Path to directory containing structure files. 130 output_dir: Optional path for per-file output. None = no file output. 131 Output files are named after the input stem, so inputs that share 132 a stem (``1crn.pdb`` and ``1crn.cif``) are rejected. Names are 133 compared without regard to ASCII case. An SDF/MOL molecule is 134 named ``stem_title`` in ``filenames`` and written under that name. 135 algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards). 136 n_points: Number of test points per atom (SR only; ignored for LR). 137 Default: 100. 138 n_slices: Number of slices per atom (LR only; ignored for SR). 139 Default: 20. 140 probe_radius: Water probe radius in Angstroms. Default: 1.4. 141 n_threads: Number of threads to use. 0 = auto-detect. Default: 0. 142 classifier: Classifier for radius assignment. None = use input radii. 143 Default: ClassifierType.CCD. With CCD, components that are not in the 144 built-in table get radii from the bond topology (``_chem_comp_bond``) 145 inside an mmCIF input, or from the bonds of an SDF/MOL molecule; there 146 is no way to pass an external CCD dictionary here (the CLI has 147 ``--ccd``). 148 include_hydrogens: Whether to include hydrogen atoms. Default: False. 149 include_hetatm: Whether to include HETATM records. Default: False. 150 151 Returns: 152 BatchDirResult with per-file details. 153 154 Raises: 155 ValueError: If input parameters are invalid, or if ``output_dir`` is 156 set and several inputs map to the same output file name. 157 FileNotFoundError: If the input directory does not exist. 158 NotADirectoryError: If ``input_dir`` is not a directory, or a parent of 159 ``output_dir`` is a file. 160 FileExistsError: If ``output_dir`` exists and is not a directory. 161 PermissionError: If the input directory cannot be read or the output 162 directory cannot be created for lack of permission. 163 OSError: If the output directory cannot be created for another reason. 164 MemoryError: If out of memory. 165 RuntimeError: For other processing errors. 166 167 Example: 168 >>> from zsasa import process_directory 169 >>> result = process_directory("path/to/pdbs/") 170 >>> print(f"Processed {result.successful}/{result.total_files} files") 171 """ 172 ffi, lib = _get_lib() 173 174 # Validate algorithm 175 if algorithm == "sr": 176 algo_int = ZSASA_ALGORITHM_SR 177 n_points_val = n_points 178 elif algorithm == "lr": 179 algo_int = ZSASA_ALGORITHM_LR 180 n_points_val = n_slices 181 else: 182 msg = f"Unknown algorithm: {algorithm}. Use 'sr' or 'lr'." 183 raise ValueError(msg) 184 185 if n_points <= 0: 186 msg = f"n_points must be positive, got {n_points}" 187 raise ValueError(msg) 188 if n_slices <= 0: 189 msg = f"n_slices must be positive, got {n_slices}" 190 raise ValueError(msg) 191 if probe_radius <= 0: 192 msg = f"probe_radius must be positive, got {probe_radius}" 193 raise ValueError(msg) 194 if n_threads < 0: 195 msg = f"n_threads must be non-negative, got {n_threads}" 196 raise ValueError(msg) 197 198 # Map classifier 199 classifier_int = int(classifier) if classifier is not None else -1 200 201 # Convert paths 202 input_dir_bytes = str(input_dir).encode("utf-8") 203 output_dir_bytes = str(output_dir).encode("utf-8") if output_dir is not None else ffi.NULL 204 205 error_code = ffi.new("int*") 206 207 handle = lib.zsasa_batch_dir_process( 208 input_dir_bytes, 209 output_dir_bytes, 210 algo_int, 211 n_points_val, 212 probe_radius, 213 n_threads, 214 classifier_int, 215 int(include_hydrogens), 216 int(include_hetatm), 217 error_code, 218 ) 219 220 if handle == ffi.NULL: 221 ec = error_code[0] 222 if ec == ZSASA_ERROR_INVALID_INPUT: 223 msg = "Invalid input parameters for directory batch processing" 224 raise ValueError(msg) 225 elif ec == ZSASA_ERROR_OUT_OF_MEMORY: 226 msg = "Out of memory during directory batch processing" 227 raise MemoryError(msg) 228 elif ec == ZSASA_ERROR_CALCULATION: 229 msg = "SASA calculation failed during directory batch processing" 230 raise RuntimeError(msg) 231 elif ec == ZSASA_ERROR_FILE_IO: 232 _raise_input_dir_error(input_dir) 233 elif ec == ZSASA_ERROR_OUTPUT_DIR: 234 _raise_output_dir_error(output_dir) 235 elif ec == ZSASA_ERROR_OUTPUT_NAME_COLLISION: 236 msg = ( 237 f"Several inputs in {input_dir} map to the same output file name " 238 "(for example the same stem with different extensions, or names that " 239 "differ only in case); split them into separate directories" 240 ) 241 raise ValueError(msg) 242 else: 243 msg = f"Directory batch processing failed with error code: {ec}" 244 raise RuntimeError(msg) 245 246 try: 247 total_files = lib.zsasa_batch_dir_get_total_files(handle) 248 successful = lib.zsasa_batch_dir_get_successful(handle) 249 failed = lib.zsasa_batch_dir_get_failed(handle) 250 251 filenames: list[str] = [] 252 n_atoms_list: list[int] = [] 253 total_sasa_list: list[float] = [] 254 status_list: list[int] = [] 255 256 for i in range(total_files): 257 fname_ptr = lib.zsasa_batch_dir_get_filename(handle, i) 258 if fname_ptr == ffi.NULL: 259 msg = ( 260 f"Internal error: zsasa_batch_dir_get_filename returned NULL " 261 f"for index {i} (total_files={total_files})" 262 ) 263 raise RuntimeError(msg) 264 filenames.append(ffi.string(fname_ptr).decode("utf-8")) 265 n_atoms_list.append(lib.zsasa_batch_dir_get_n_atoms(handle, i)) 266 sasa = lib.zsasa_batch_dir_get_total_sasa(handle, i) 267 total_sasa_list.append(float("nan") if math.isnan(sasa) else sasa) 268 st = lib.zsasa_batch_dir_get_status(handle, i) 269 if st not in (0, 1): 270 msg = f"Internal error: unexpected status {st} for file index {i} (expected 0 or 1)" 271 raise RuntimeError(msg) 272 status_list.append(st) 273 274 return BatchDirResult( 275 total_files=total_files, 276 successful=successful, 277 failed=failed, 278 filenames=filenames, 279 n_atoms=n_atoms_list, 280 total_sasa=total_sasa_list, 281 status=status_list, 282 ) 283 finally: 284 lib.zsasa_batch_dir_free(handle)
Process all supported structure files in a directory for SASA calculation.
Supported formats: PDB (.pdb, .ent), mmCIF (.cif, .mmcif), BinaryCIF (.bcif), JSON (.json), SDF/MOL (.sdf, .mol), and their gzip- or zstd-compressed variants (.gz, .zst).
Args:
input_dir: Path to directory containing structure files.
output_dir: Optional path for per-file output. None = no file output.
Output files are named after the input stem, so inputs that share
a stem (1crn.pdb and 1crn.cif) are rejected. Names are
compared without regard to ASCII case. An SDF/MOL molecule is
named stem_title in filenames and written under that name.
algorithm: Algorithm to use: "sr" (Shrake-Rupley) or "lr" (Lee-Richards).
n_points: Number of test points per atom (SR only; ignored for LR).
Default: 100.
n_slices: Number of slices per atom (LR only; ignored for SR).
Default: 20.
probe_radius: Water probe radius in Angstroms. Default: 1.4.
n_threads: Number of threads to use. 0 = auto-detect. Default: 0.
classifier: Classifier for radius assignment. None = use input radii.
Default: ClassifierType.CCD. With CCD, components that are not in the
built-in table get radii from the bond topology (_chem_comp_bond)
inside an mmCIF input, or from the bonds of an SDF/MOL molecule; there
is no way to pass an external CCD dictionary here (the CLI has
--ccd).
include_hydrogens: Whether to include hydrogen atoms. Default: False.
include_hetatm: Whether to include HETATM records. Default: False.
Returns: BatchDirResult with per-file details.
Raises:
ValueError: If input parameters are invalid, or if output_dir is
set and several inputs map to the same output file name.
FileNotFoundError: If the input directory does not exist.
NotADirectoryError: If input_dir is not a directory, or a parent of
output_dir is a file.
FileExistsError: If output_dir exists and is not a directory.
PermissionError: If the input directory cannot be read or the output
directory cannot be created for lack of permission.
OSError: If the output directory cannot be created for another reason.
MemoryError: If out of memory.
RuntimeError: For other processing errors.
Example:
from zsasa import process_directory result = process_directory("path/to/pdbs/") print(f"Processed {result.successful}/{result.total_files} files")
61@dataclass 62class BatchDirResult: 63 """Result of directory batch processing. 64 65 Attributes: 66 total_files: Number of supported structure files found in the directory. 67 successful: Number of successfully processed files. 68 failed: Number of files that failed processing. 69 filenames: List of filenames (file name only, without directory path). 70 n_atoms: Per-file atom counts. 71 total_sasa: Per-file total SASA in Angstroms² (NaN for failed files). 72 status: Per-file status (1=ok, 0=failed). 73 """ 74 75 total_files: int 76 successful: int 77 failed: int 78 filenames: list[str] 79 n_atoms: list[int] 80 total_sasa: list[float] 81 status: list[int] 82 83 def __post_init__(self) -> None: 84 n = len(self.filenames) 85 if len(self.n_atoms) != n or len(self.total_sasa) != n or len(self.status) != n: 86 msg = ( 87 f"All per-file lists must have the same length as filenames ({n}), " 88 f"got n_atoms={len(self.n_atoms)}, total_sasa={len(self.total_sasa)}, " 89 f"status={len(self.status)}" 90 ) 91 raise ValueError(msg) 92 if self.total_files != n: 93 msg = f"total_files ({self.total_files}) != len(filenames) ({n})" 94 raise ValueError(msg) 95 if self.successful + self.failed != self.total_files: 96 msg = ( 97 f"successful ({self.successful}) + failed ({self.failed}) " 98 f"!= total_files ({self.total_files})" 99 ) 100 raise ValueError(msg) 101 102 def __repr__(self) -> str: 103 return ( 104 f"BatchDirResult(total_files={self.total_files}, " 105 f"successful={self.successful}, failed={self.failed})" 106 )
Result of directory batch processing.
Attributes: total_files: Number of supported structure files found in the directory. successful: Number of successfully processed files. failed: Number of files that failed processing. filenames: List of filenames (file name only, without directory path). n_atoms: Per-file atom counts. total_sasa: Per-file total SASA in Angstroms² (NaN for failed files). status: Per-file status (1=ok, 0=failed).
26class ClassifierType(IntEnum): 27 """Available classifier types for atom radius assignment. 28 29 Attributes: 30 CCD: The default classifier. The CLI and ``process_directory`` derive radii 31 from bond topology for components that are not in the built-in table 32 (inline ``_chem_comp_bond`` data of an mmCIF file, SDF bonds, and with 33 the CLI an external ``--ccd`` dictionary). Everything that goes through 34 the C classify API does not: ``get_radius``, ``get_atom_class``, 35 ``classify_atoms`` and the gemmi, BioPython and Biotite integrations use 36 the built-in ProtOr table only, with no bond-topology analysis of 37 unknown components. Such atoms get a NaN radius (``None`` from 38 ``get_radius``); the integrations fill them from the element. 39 PROTOR: Static ProtOr-compatible radii without CCD resource parsing. A 40 separate value, not an alias of CCD; for ``get_radius`` and 41 ``classify_atoms`` it returns the same radii as CCD. 42 NACCESS: NACCESS-compatible radii. 43 OONS: Ooi, Oobatake, Nemethy, Scheraga radii. 44 """ 45 46 NACCESS = ZSASA_CLASSIFIER_NACCESS 47 PROTOR = ZSASA_CLASSIFIER_PROTOR 48 OONS = ZSASA_CLASSIFIER_OONS 49 CCD = ZSASA_CLASSIFIER_CCD
Available classifier types for atom radius assignment.
Attributes:
CCD: The default classifier. The CLI and process_directory derive radii
from bond topology for components that are not in the built-in table
(inline _chem_comp_bond data of an mmCIF file, SDF bonds, and with
the CLI an external --ccd dictionary). Everything that goes through
the C classify API does not: get_radius, get_atom_class,
classify_atoms and the gemmi, BioPython and Biotite integrations use
the built-in ProtOr table only, with no bond-topology analysis of
unknown components. Such atoms get a NaN radius (None from
get_radius); the integrations fill them from the element.
PROTOR: Static ProtOr-compatible radii without CCD resource parsing. A
separate value, not an alias of CCD; for get_radius and
classify_atoms it returns the same radii as CCD.
NACCESS: NACCESS-compatible radii.
OONS: Ooi, Oobatake, Nemethy, Scheraga radii.
52class AtomClass(IntEnum): 53 """Atom polarity classes. 54 55 Attributes: 56 POLAR: Polar atoms (e.g., N, O). 57 APOLAR: Apolar/hydrophobic atoms (e.g., C). 58 UNKNOWN: Unknown classification. 59 """ 60 61 POLAR = ZSASA_ATOM_CLASS_POLAR 62 APOLAR = ZSASA_ATOM_CLASS_APOLAR 63 UNKNOWN = ZSASA_ATOM_CLASS_UNKNOWN
Atom polarity classes.
Attributes: POLAR: Polar atoms (e.g., N, O). APOLAR: Apolar/hydrophobic atoms (e.g., C). UNKNOWN: Unknown classification.
203@dataclass 204class ClassificationResult: 205 """Result of batch atom classification. 206 207 Attributes: 208 radii: Per-atom van der Waals radii in Angstroms. 209 NaN values indicate atoms not found in the classifier. 210 Use np.isnan(result.radii) to find unknown atoms. 211 classes: Per-atom polarity classes (AtomClass constants). 212 213 Example: 214 >>> result = classify_atoms(residues, atoms) 215 >>> unknown_mask = np.isnan(result.radii) 216 >>> if unknown_mask.any(): 217 ... print(f"Found {unknown_mask.sum()} unknown atoms") 218 """ 219 220 radii: NDArray[np.float64] 221 classes: NDArray[np.int32] 222 223 def __repr__(self) -> str: 224 return f"ClassificationResult(n_atoms={len(self.radii)})"
Result of batch atom classification.
Attributes: radii: Per-atom van der Waals radii in Angstroms. NaN values indicate atoms not found in the classifier. Use np.isnan(result.radii) to find unknown atoms. classes: Per-atom polarity classes (AtomClass constants).
Example:
result = classify_atoms(residues, atoms) unknown_mask = np.isnan(result.radii) if unknown_mask.any(): ... print(f"Found {unknown_mask.sum()} unknown atoms")
66def get_radius( 67 residue: str, 68 atom: str, 69 classifier_type: ClassifierType = ClassifierType.CCD, 70) -> float | None: 71 """Get van der Waals radius for an atom using the specified classifier. 72 73 Args: 74 residue: Residue name (e.g., "ALA", "GLY"). 75 atom: Atom name (e.g., "CA", "CB"). 76 classifier_type: Classifier to use. Default: ClassifierType.CCD. 77 78 Returns: 79 Radius in Angstroms, or None if atom is not found in classifier. 80 81 Note: 82 ``ClassifierType.CCD`` looks the atom up in the built-in ProtOr table only. 83 This function cannot load chemical component definitions, so a ligand or 84 modified residue that is not in the table gives None, where the CLI would 85 analyze its bond topology. 86 87 NACCESS and OONS have no entries for hydrogens or ligands. For those 88 atoms the element is guessed from the names: a name starting with 89 H, C, N, O, P or S is that element ("HG" is hydrogen, "NA" in "HEM" 90 is nitrogen), and an ion is recognized by a residue name equal to its 91 atom name ("ZN" in "ZN"). A metal named like an organic atom inside a 92 larger residue cannot be recognized this way, so prefer 93 ``guess_radius(element)`` when the element is known. 94 95 Example: 96 >>> from zsasa import get_radius, ClassifierType 97 >>> get_radius("ALA", "CA") 98 1.88 99 >>> get_radius("ALA", "CA", ClassifierType.NACCESS) 100 1.87 101 >>> get_radius("ALA", "XX") # Unknown atom 102 None 103 """ 104 _, lib = _get_lib() 105 radius = lib.zsasa_classifier_get_radius( 106 classifier_type, 107 residue.encode("utf-8"), 108 atom.encode("utf-8"), 109 ) 110 if np.isnan(radius): 111 return None 112 return radius
Get van der Waals radius for an atom using the specified classifier.
Args: residue: Residue name (e.g., "ALA", "GLY"). atom: Atom name (e.g., "CA", "CB"). classifier_type: Classifier to use. Default: ClassifierType.CCD.
Returns: Radius in Angstroms, or None if atom is not found in classifier.
Note:
ClassifierType.CCD looks the atom up in the built-in ProtOr table only.
This function cannot load chemical component definitions, so a ligand or
modified residue that is not in the table gives None, where the CLI would
analyze its bond topology.
NACCESS and OONS have no entries for hydrogens or ligands. For those
atoms the element is guessed from the names: a name starting with
H, C, N, O, P or S is that element ("HG" is hydrogen, "NA" in "HEM"
is nitrogen), and an ion is recognized by a residue name equal to its
atom name ("ZN" in "ZN"). A metal named like an organic atom inside a
larger residue cannot be recognized this way, so prefer
``guess_radius(element)`` when the element is known.
Example:
from zsasa import get_radius, ClassifierType get_radius("ALA", "CA") 1.88 get_radius("ALA", "CA", ClassifierType.NACCESS) 1.87 get_radius("ALA", "XX") # Unknown atom None
115def get_atom_class( 116 residue: str, 117 atom: str, 118 classifier_type: ClassifierType = ClassifierType.CCD, 119) -> int: 120 """Get atom polarity class using the specified classifier. 121 122 Args: 123 residue: Residue name (e.g., "ALA", "GLY"). 124 atom: Atom name (e.g., "CA", "CB"). 125 classifier_type: Classifier to use. Default: ClassifierType.CCD. 126 127 Returns: 128 AtomClass constant (POLAR, APOLAR, or UNKNOWN). ``ClassifierType.CCD`` 129 uses the built-in ProtOr table only (see ``classify_atoms``), so atoms of 130 components outside it are UNKNOWN. 131 132 Example: 133 >>> from zsasa import get_atom_class, AtomClass 134 >>> get_atom_class("ALA", "CA") == AtomClass.APOLAR 135 True 136 >>> get_atom_class("ALA", "O") == AtomClass.POLAR 137 True 138 """ 139 _, lib = _get_lib() 140 return lib.zsasa_classifier_get_class( 141 classifier_type, 142 residue.encode("utf-8"), 143 atom.encode("utf-8"), 144 )
Get atom polarity class using the specified classifier.
Args: residue: Residue name (e.g., "ALA", "GLY"). atom: Atom name (e.g., "CA", "CB"). classifier_type: Classifier to use. Default: ClassifierType.CCD.
Returns:
AtomClass constant (POLAR, APOLAR, or UNKNOWN). ClassifierType.CCD
uses the built-in ProtOr table only (see classify_atoms), so atoms of
components outside it are UNKNOWN.
Example:
from zsasa import get_atom_class, AtomClass get_atom_class("ALA", "CA") == AtomClass.APOLAR True get_atom_class("ALA", "O") == AtomClass.POLAR True
147def guess_radius(element: str) -> float | None: 148 """Guess van der Waals radius from element symbol. 149 150 Args: 151 element: Element symbol (e.g., "C", "N", "FE"). 152 Case-insensitive, whitespace is trimmed. 153 154 Returns: 155 Radius in Angstroms, or None if element is not recognized. 156 157 Example: 158 >>> from zsasa import guess_radius 159 >>> guess_radius("C") 160 1.7 161 >>> guess_radius("FE") 162 1.26 163 >>> guess_radius("XX") # Unknown 164 None 165 """ 166 _, lib = _get_lib() 167 radius = lib.zsasa_guess_radius(element.encode("utf-8")) 168 if np.isnan(radius): 169 return None 170 return radius
Guess van der Waals radius from element symbol.
Args: element: Element symbol (e.g., "C", "N", "FE"). Case-insensitive, whitespace is trimmed.
Returns: Radius in Angstroms, or None if element is not recognized.
Example:
from zsasa import guess_radius guess_radius("C") 1.7 guess_radius("FE") 1.26 guess_radius("XX") # Unknown None
173def guess_radius_from_atom_name(atom_name: str) -> float | None: 174 """Guess van der Waals radius from PDB-style atom name. 175 176 Extracts element symbol from atom name following PDB conventions: 177 - Leading space indicates single-char element (e.g., " CA " = Carbon alpha) 178 - No leading space may indicate 2-char element (e.g., "FE " = Iron) 179 180 Pass the name with its column padding: a trimmed "CA" or "HG" is read as 181 calcium or mercury. 182 183 Args: 184 atom_name: PDB-style atom name (e.g., " CA ", "FE "). 185 186 Returns: 187 Radius in Angstroms, or None if element cannot be determined. 188 189 Example: 190 >>> from zsasa import guess_radius_from_atom_name 191 >>> guess_radius_from_atom_name(" CA ") # Carbon alpha 192 1.7 193 >>> guess_radius_from_atom_name("FE ") # Iron 194 1.26 195 """ 196 _, lib = _get_lib() 197 radius = lib.zsasa_guess_radius_from_atom_name(atom_name.encode("utf-8")) 198 if np.isnan(radius): 199 return None 200 return radius
Guess van der Waals radius from PDB-style atom name.
Extracts element symbol from atom name following PDB conventions:
- Leading space indicates single-char element (e.g., " CA " = Carbon alpha)
- No leading space may indicate 2-char element (e.g., "FE " = Iron)
Pass the name with its column padding: a trimmed "CA" or "HG" is read as calcium or mercury.
Args: atom_name: PDB-style atom name (e.g., " CA ", "FE ").
Returns: Radius in Angstroms, or None if element cannot be determined.
Example:
from zsasa import guess_radius_from_atom_name guess_radius_from_atom_name(" CA ") # Carbon alpha 1.7 guess_radius_from_atom_name("FE ") # Iron 1.26
227def classify_atoms( 228 residues: list[str], 229 atoms: list[str], 230 classifier_type: ClassifierType = ClassifierType.CCD, 231 *, 232 include_classes: bool = True, 233) -> ClassificationResult: 234 """Classify multiple atoms at once (batch operation). 235 236 This is more efficient than calling get_radius for each atom individually. 237 238 ``ClassifierType.CCD`` (the default) and ``ClassifierType.PROTOR`` both look the 239 atoms up in the built-in ProtOr table. There is no bond-topology analysis of 240 components that are not in the table: the C classify API has no source for 241 chemical component definitions, so those atoms get a NaN radius and UNKNOWN 242 class. The CLI and ``process_directory`` do analyze them (see ``ClassifierType``). 243 244 Args: 245 residues: List of residue names. 246 atoms: List of atom names (must be same length as residues). 247 classifier_type: Classifier to use. Default: ClassifierType.CCD. 248 include_classes: Whether to compute atom classes. Default: True. 249 250 Returns: 251 ClassificationResult with radii and classes arrays. 252 Unknown atoms have NaN radius and UNKNOWN class. With NACCESS and 253 OONS, atoms outside the tables have UNKNOWN class and a radius 254 guessed from the names where possible (see get_radius). 255 256 Raises: 257 ValueError: If residues and atoms have different lengths. 258 259 Example: 260 >>> from zsasa import classify_atoms 261 >>> result = classify_atoms( 262 ... ["ALA", "ALA", "GLY"], 263 ... ["CA", "O", "N"], 264 ... ) 265 >>> result.radii 266 array([1.88, 1.42, 1.64]) 267 """ 268 ffi, lib = _get_lib() 269 270 if len(residues) != len(atoms): 271 msg = f"residues and atoms must have same length: {len(residues)} != {len(atoms)}" 272 raise ValueError(msg) 273 274 n_atoms = len(residues) 275 if n_atoms == 0: 276 return ClassificationResult( 277 radii=np.array([], dtype=np.float64), 278 classes=np.array([], dtype=np.int32), 279 ) 280 281 # Encode strings and create cffi arrays 282 # Keep references to prevent garbage collection 283 residues_bytes = [ffi.new("char[]", r.encode("utf-8")) for r in residues] 284 atoms_bytes = [ffi.new("char[]", a.encode("utf-8")) for a in atoms] 285 286 residues_arr = ffi.new("char*[]", residues_bytes) 287 atoms_arr = ffi.new("char*[]", atoms_bytes) 288 289 # Allocate output arrays 290 radii = np.zeros(n_atoms, dtype=np.float64) 291 classes = np.zeros(n_atoms, dtype=np.int32) if include_classes else None 292 293 radii_ptr = ffi.cast("double*", radii.ctypes.data) 294 classes_ptr = ffi.cast("int*", classes.ctypes.data) if include_classes else ffi.NULL 295 296 result = lib.zsasa_classify_atoms( 297 classifier_type, 298 residues_arr, 299 atoms_arr, 300 n_atoms, 301 radii_ptr, 302 classes_ptr, 303 ) 304 305 if result == ZSASA_ERROR_INVALID_INPUT: 306 msg = f"Invalid classifier type: {classifier_type}" 307 raise ValueError(msg) 308 elif result != ZSASA_OK: 309 msg = f"Classification error: {result}" 310 raise RuntimeError(msg) 311 312 if not include_classes: 313 classes = np.full(n_atoms, AtomClass.UNKNOWN, dtype=np.int32) 314 315 return ClassificationResult(radii=radii, classes=classes)
Classify multiple atoms at once (batch operation).
This is more efficient than calling get_radius for each atom individually.
ClassifierType.CCD (the default) and ClassifierType.PROTOR both look the
atoms up in the built-in ProtOr table. There is no bond-topology analysis of
components that are not in the table: the C classify API has no source for
chemical component definitions, so those atoms get a NaN radius and UNKNOWN
class. The CLI and process_directory do analyze them (see ClassifierType).
Args: residues: List of residue names. atoms: List of atom names (must be same length as residues). classifier_type: Classifier to use. Default: ClassifierType.CCD. include_classes: Whether to compute atom classes. Default: True.
Returns: ClassificationResult with radii and classes arrays. Unknown atoms have NaN radius and UNKNOWN class. With NACCESS and OONS, atoms outside the tables have UNKNOWN class and a radius guessed from the names where possible (see get_radius).
Raises: ValueError: If residues and atoms have different lengths.
Example:
from zsasa import classify_atoms result = classify_atoms( ... ["ALA", "ALA", "GLY"], ... ["CA", "O", "N"], ... ) result.radii array([1.88, 1.42, 1.64])
36def get_max_sasa(residue_name: str) -> float | None: 37 """Get maximum SASA value for a standard amino acid. 38 39 Values from Tien et al. (2013) "Maximum allowed solvent accessibilities 40 of residues in proteins". 41 42 Args: 43 residue_name: 3-letter residue code (e.g., "ALA", "GLY"). 44 45 Returns: 46 Maximum SASA in Angstroms², or None if residue is not a standard amino acid. 47 48 Example: 49 >>> from zsasa import get_max_sasa 50 >>> get_max_sasa("ALA") 51 129.0 52 >>> get_max_sasa("TRP") 53 285.0 54 >>> get_max_sasa("HOH") # Water - not a standard amino acid 55 None 56 """ 57 _, lib = _get_lib() 58 max_sasa = lib.zsasa_get_max_sasa(residue_name.encode("utf-8")) 59 if np.isnan(max_sasa): 60 return None 61 return max_sasa
Get maximum SASA value for a standard amino acid.
Values from Tien et al. (2013) "Maximum allowed solvent accessibilities of residues in proteins".
Args: residue_name: 3-letter residue code (e.g., "ALA", "GLY").
Returns: Maximum SASA in Angstroms², or None if residue is not a standard amino acid.
Example:
from zsasa import get_max_sasa get_max_sasa("ALA") 129.0 get_max_sasa("TRP") 285.0 get_max_sasa("HOH") # Water - not a standard amino acid None
64def calculate_rsa(sasa: float, residue_name: str) -> float | None: 65 """Calculate RSA (Relative Solvent Accessibility) for a single residue. 66 67 RSA = SASA / MaxSASA 68 69 Args: 70 sasa: Observed SASA value in Angstroms². 71 residue_name: 3-letter residue code (e.g., "ALA", "GLY"). 72 73 Returns: 74 RSA value (typically 0.0-1.0), or None if residue is not a standard amino acid. 75 Note: RSA > 1.0 is possible for exposed terminal residues. 76 77 Example: 78 >>> from zsasa import calculate_rsa 79 >>> calculate_rsa(64.5, "ALA") # 64.5 / 129.0 = 0.5 80 0.5 81 >>> calculate_rsa(150.0, "GLY") # RSA > 1.0 is possible 82 1.4423076923076923 83 """ 84 _, lib = _get_lib() 85 rsa = lib.zsasa_calculate_rsa(sasa, residue_name.encode("utf-8")) 86 if np.isnan(rsa): 87 return None 88 return rsa
Calculate RSA (Relative Solvent Accessibility) for a single residue.
RSA = SASA / MaxSASA
Args: sasa: Observed SASA value in Angstroms². residue_name: 3-letter residue code (e.g., "ALA", "GLY").
Returns: RSA value (typically 0.0-1.0), or None if residue is not a standard amino acid. Note: RSA > 1.0 is possible for exposed terminal residues.
Example:
from zsasa import calculate_rsa calculate_rsa(64.5, "ALA") # 64.5 / 129.0 = 0.5 0.5 calculate_rsa(150.0, "GLY") # RSA > 1.0 is possible 1.4423076923076923
91def calculate_rsa_batch( 92 sasas: NDArray[np.float64] | list[float], 93 residue_names: list[str], 94) -> NDArray[np.float64]: 95 """Calculate RSA for multiple residues at once (batch operation). 96 97 This is more efficient than calling calculate_rsa for each residue individually. 98 99 Args: 100 sasas: Array of SASA values in Angstroms². 101 residue_names: List of 3-letter residue codes (must be same length as sasas). 102 103 Returns: 104 Array of RSA values. NaN values indicate non-standard amino acids. 105 106 Raises: 107 ValueError: If sasas and residue_names have different lengths. 108 109 Example: 110 >>> import numpy as np 111 >>> from zsasa import calculate_rsa_batch 112 >>> sasas = np.array([64.5, 52.0, 100.0]) 113 >>> residues = ["ALA", "GLY", "HOH"] # HOH is not standard 114 >>> rsa = calculate_rsa_batch(sasas, residues) 115 >>> rsa 116 array([0.5 , 0.5 , nan]) 117 """ 118 ffi, lib = _get_lib() 119 120 sasas = np.ascontiguousarray(sasas, dtype=np.float64) 121 n_residues = len(sasas) 122 123 if len(residue_names) != n_residues: 124 msg = f"sasas and residue_names must have same length: {n_residues} != {len(residue_names)}" 125 raise ValueError(msg) 126 127 if n_residues == 0: 128 return np.array([], dtype=np.float64) 129 130 # Encode strings and create cffi array 131 # Keep references to prevent garbage collection 132 residues_bytes = [ffi.new("char[]", r.encode("utf-8")) for r in residue_names] 133 residues_arr = ffi.new("char*[]", residues_bytes) 134 135 # Allocate output array 136 rsa_out = np.zeros(n_residues, dtype=np.float64) 137 138 sasas_ptr = ffi.cast("double*", sasas.ctypes.data) 139 rsa_ptr = ffi.cast("double*", rsa_out.ctypes.data) 140 141 result = lib.zsasa_calculate_rsa_batch(sasas_ptr, residues_arr, n_residues, rsa_ptr) 142 if result != ZSASA_OK: 143 msg = f"RSA batch calculation failed with error code: {result}" 144 raise RuntimeError(msg) 145 146 return rsa_out
Calculate RSA for multiple residues at once (batch operation).
This is more efficient than calling calculate_rsa for each residue individually.
Args: sasas: Array of SASA values in Angstroms². residue_names: List of 3-letter residue codes (must be same length as sasas).
Returns: Array of RSA values. NaN values indicate non-standard amino acids.
Raises: ValueError: If sasas and residue_names have different lengths.
Example:
import numpy as np from zsasa import calculate_rsa_batch sasas = np.array([64.5, 52.0, 100.0]) residues = ["ALA", "GLY", "HOH"] # HOH is not standard rsa = calculate_rsa_batch(sasas, residues) rsa array([0.5 , 0.5 , nan])
41@dataclass 42class ResidueResult: 43 """Per-residue SASA result. 44 45 Attributes: 46 chain_id: Chain identifier (e.g., "A", "B"). 47 residue_id: Residue sequence number. 48 residue_name: 3-letter residue name (e.g., "ALA", "GLY"). 49 total_area: Total SASA for this residue in A^2. 50 polar_area: Polar SASA (N, O atoms, etc.) in A^2. 51 apolar_area: Apolar SASA (C atoms, etc.) in A^2. 52 rsa: Relative Solvent Accessibility (0.0-1.0+), or None for 53 non-standard amino acids. 54 n_atoms: Number of atoms in this residue. 55 insertion_code: Insertion code of the residue, ``""`` for a residue 56 without one. Residues ``10``, ``10A`` and ``10B`` have the same 57 ``residue_id`` and differ only here. 58 """ 59 60 chain_id: str 61 residue_id: int 62 residue_name: str 63 total_area: float 64 polar_area: float 65 apolar_area: float 66 rsa: float | None 67 n_atoms: int 68 insertion_code: str = "" 69 70 def __repr__(self) -> str: 71 rsa_str = f"{self.rsa:.3f}" if self.rsa is not None else "None" 72 residue = f"{self.residue_name}{self.residue_id}{self.insertion_code}" 73 return ( 74 f"ResidueResult({self.chain_id}:{residue}, " 75 f"total={self.total_area:.1f}, rsa={rsa_str}, n_atoms={self.n_atoms})" 76 )
Per-residue SASA result.
Attributes:
chain_id: Chain identifier (e.g., "A", "B").
residue_id: Residue sequence number.
residue_name: 3-letter residue name (e.g., "ALA", "GLY").
total_area: Total SASA for this residue in A^2.
polar_area: Polar SASA (N, O atoms, etc.) in A^2.
apolar_area: Apolar SASA (C atoms, etc.) in A^2.
rsa: Relative Solvent Accessibility (0.0-1.0+), or None for
non-standard amino acids.
n_atoms: Number of atoms in this residue.
insertion_code: Insertion code of the residue, "" for a residue
without one. Residues 10, 10A and 10B have the same
residue_id and differ only here.
79def aggregate_by_residue( 80 atom_areas: NDArray[np.float64], 81 chain_ids: list[str], 82 residue_ids: list[int], 83 residue_names: list[str], 84 atom_classes: NDArray[np.int32] | None = None, 85 insertion_codes: list[str] | None = None, 86) -> list[ResidueResult]: 87 """Aggregate per-atom SASA values to per-residue. 88 89 Groups atoms by (chain_id, residue_id, insertion_code) and sums their SASA 90 values. Also calculates polar/apolar breakdown and RSA if atom classes are 91 provided. 92 93 Pass ``insertion_codes`` for structures that use them (antibody numbering, 94 for example): without them, residues such as ``10``, ``10A`` and ``10B`` 95 have the same ``residue_id`` and are merged into one result. 96 97 Args: 98 atom_areas: Per-atom SASA values in A^2. 99 chain_ids: Chain ID for each atom. 100 residue_ids: Residue sequence number for each atom. 101 residue_names: Residue name for each atom. 102 atom_classes: Optional per-atom polarity classes (AtomClass values). 103 If provided, polar_area and apolar_area will be calculated. 104 insertion_codes: Optional insertion code of the residue of each atom, 105 ``""`` (or a blank) for a residue without one. If 106 omitted, no residue has an insertion code. 107 108 Returns: 109 List of ResidueResult objects, one per unique residue. 110 Results are ordered by appearance in the input (preserves chain/residue order). 111 112 Example: 113 >>> import numpy as np 114 >>> from zsasa.analysis import aggregate_by_residue 115 >>> 116 >>> atom_areas = np.array([10.0, 20.0, 15.0, 25.0]) 117 >>> chain_ids = ["A", "A", "A", "A"] 118 >>> residue_ids = [1, 1, 2, 2] 119 >>> residue_names = ["ALA", "ALA", "GLY", "GLY"] 120 >>> 121 >>> residues = aggregate_by_residue( 122 ... atom_areas, chain_ids, residue_ids, residue_names 123 ... ) 124 >>> len(residues) 125 2 126 >>> residues[0].total_area 127 30.0 128 """ 129 n_atoms = len(atom_areas) 130 if n_atoms == 0: 131 return [] 132 133 # Validate input lengths 134 if len(chain_ids) != n_atoms: 135 msg = f"chain_ids length ({len(chain_ids)}) != atom_areas length ({n_atoms})" 136 raise ValueError(msg) 137 if len(residue_ids) != n_atoms: 138 msg = f"residue_ids length ({len(residue_ids)}) != atom_areas length ({n_atoms})" 139 raise ValueError(msg) 140 if len(residue_names) != n_atoms: 141 msg = f"residue_names length ({len(residue_names)}) != atom_areas length ({n_atoms})" 142 raise ValueError(msg) 143 if atom_classes is not None and len(atom_classes) != n_atoms: 144 msg = f"atom_classes length ({len(atom_classes)}) != atom_areas length ({n_atoms})" 145 raise ValueError(msg) 146 if insertion_codes is not None and len(insertion_codes) != n_atoms: 147 msg = f"insertion_codes length ({len(insertion_codes)}) != atom_areas length ({n_atoms})" 148 raise ValueError(msg) 149 150 # Group atoms by (chain_id, residue_id, insertion_code) 151 # Use dict to preserve insertion order (Python 3.7+) 152 residue_data: dict[tuple[str, int, str], dict] = {} 153 154 for i in range(n_atoms): 155 # Structure libraries write a blank for "no insertion code" 156 insertion_code = insertion_codes[i].strip() if insertion_codes is not None else "" 157 key = (chain_ids[i], residue_ids[i], insertion_code) 158 159 if key not in residue_data: 160 residue_data[key] = { 161 "residue_name": residue_names[i], 162 "total_area": 0.0, 163 "polar_area": 0.0, 164 "apolar_area": 0.0, 165 "n_atoms": 0, 166 } 167 168 data = residue_data[key] 169 area = float(atom_areas[i]) 170 data["total_area"] += area 171 data["n_atoms"] += 1 172 173 if atom_classes is not None: 174 atom_class = atom_classes[i] 175 if atom_class == AtomClass.POLAR: 176 data["polar_area"] += area 177 elif atom_class == AtomClass.APOLAR: 178 data["apolar_area"] += area 179 180 # Build result list 181 results = [] 182 for (chain_id, residue_id, insertion_code), data in residue_data.items(): 183 residue_name = data["residue_name"] 184 total_area = data["total_area"] 185 186 # Calculate RSA if this is a standard amino acid 187 max_sasa = MAX_SASA.get(residue_name) 188 rsa = total_area / max_sasa if max_sasa is not None else None 189 190 results.append( 191 ResidueResult( 192 chain_id=chain_id, 193 residue_id=residue_id, 194 residue_name=residue_name, 195 total_area=total_area, 196 polar_area=data["polar_area"], 197 apolar_area=data["apolar_area"], 198 rsa=rsa, 199 n_atoms=data["n_atoms"], 200 insertion_code=insertion_code, 201 ) 202 ) 203 204 return results
Aggregate per-atom SASA values to per-residue.
Groups atoms by (chain_id, residue_id, insertion_code) and sums their SASA values. Also calculates polar/apolar breakdown and RSA if atom classes are provided.
Pass insertion_codes for structures that use them (antibody numbering,
for example): without them, residues such as 10, 10A and 10B
have the same residue_id and are merged into one result.
Args:
atom_areas: Per-atom SASA values in A^2.
chain_ids: Chain ID for each atom.
residue_ids: Residue sequence number for each atom.
residue_names: Residue name for each atom.
atom_classes: Optional per-atom polarity classes (AtomClass values).
If provided, polar_area and apolar_area will be calculated.
insertion_codes: Optional insertion code of the residue of each atom,
"" (or a blank) for a residue without one. If
omitted, no residue has an insertion code.
Returns: List of ResidueResult objects, one per unique residue. Results are ordered by appearance in the input (preserves chain/residue order).
Example:
import numpy as np from zsasa.analysis import aggregate_by_residue
atom_areas = np.array([10.0, 20.0, 15.0, 25.0]) chain_ids = ["A", "A", "A", "A"] residue_ids = [1, 1, 2, 2] residue_names = ["ALA", "ALA", "GLY", "GLY"]
residues = aggregate_by_residue( ... atom_areas, chain_ids, residue_ids, residue_names ... ) len(residues) 2 residues[0].total_area 30.0
207def aggregate_from_result(result: SasaResultWithAtoms) -> list[ResidueResult]: 208 """Aggregate per-atom SASA from a SasaResultWithAtoms to per-residue. 209 210 This is a convenience wrapper that extracts the necessary data from 211 a SasaResultWithAtoms object returned by gemmi integration functions. 212 213 Args: 214 result: A SasaResultWithAtoms object from calculate_sasa_from_structure 215 or calculate_sasa_from_model. 216 217 Returns: 218 List of ResidueResult objects, one per unique residue. 219 220 Example: 221 >>> from zsasa.integrations.gemmi import calculate_sasa_from_structure 222 >>> from zsasa.analysis import aggregate_from_result 223 >>> 224 >>> result = calculate_sasa_from_structure("protein.cif") 225 >>> residues = aggregate_from_result(result) 226 >>> 227 >>> # Print buried residues (RSA < 0.25) 228 >>> for res in residues: 229 ... if res.rsa is not None and res.rsa < 0.25: 230 ... print(f"{res.chain_id}:{res.residue_name}{res.residue_id}: {res.rsa:.1%}") 231 """ 232 return aggregate_by_residue( 233 atom_areas=result.atom_areas, 234 chain_ids=result.atom_data.chain_ids, 235 residue_ids=result.atom_data.residue_ids, 236 residue_names=result.atom_data.residue_names, 237 atom_classes=result.atom_classes, 238 insertion_codes=result.atom_data.insertion_codes, 239 )
Aggregate per-atom SASA from a SasaResultWithAtoms to per-residue.
This is a convenience wrapper that extracts the necessary data from a SasaResultWithAtoms object returned by gemmi integration functions.
Args: result: A SasaResultWithAtoms object from calculate_sasa_from_structure or calculate_sasa_from_model.
Returns: List of ResidueResult objects, one per unique residue.
Example:
from zsasa.integrations.gemmi import calculate_sasa_from_structure from zsasa.analysis import aggregate_from_result
result = calculate_sasa_from_structure("protein.cif") residues = aggregate_from_result(result)
Print buried residues (RSA < 0.25)
for res in residues: ... if res.rsa is not None and res.rsa < 0.25: ... print(f"{res.chain_id}:{res.residue_name}{res.residue_id}: {res.rsa:.1%}")
393def get_version() -> str: 394 """Get the library version string.""" 395 ffi, lib = _get_lib() 396 return ffi.string(lib.zsasa_version()).decode("utf-8")
Get the library version string.