0
# Coordinate Transformations
1
2
MDAnalysis provides comprehensive tools for geometric transformations and coordinate manipulations. These capabilities enable structural alignment, coordinate system transformations, and periodic boundary condition handling essential for molecular dynamics analysis.
3
4
## Overview
5
6
Coordinate transformations in MDAnalysis fall into several categories:
7
8
- **Basic Transformations**: Translation, rotation, scaling
9
- **Structural Alignment**: Fitting molecules to reference structures
10
- **Periodic Boundary Handling**: Wrapping and unwrapping coordinates
11
- **On-the-fly Transformations**: Trajectory transformations during reading
12
- **Reference Fitting**: Removing overall motion and rotation
13
14
## Basic Transformations
15
16
### Translation
17
18
```python { .api }
19
def translate(atomgroup, t):
20
"""
21
Translate atomic coordinates by a vector.
22
23
Parameters
24
----------
25
atomgroup : AtomGroup
26
Atoms to translate.
27
t : array-like
28
Translation vector as 3-element array [dx, dy, dz] in Angstrom.
29
30
Examples
31
--------
32
>>> # Move protein 10 Å along x-axis
33
>>> protein.translate([10.0, 0.0, 0.0])
34
35
>>> # Center protein at origin
36
>>> com = protein.center_of_mass()
37
>>> protein.translate(-com)
38
39
>>> # Move to specific position
40
>>> target_position = [50.0, 25.0, 75.0]
41
>>> current_com = protein.center_of_mass()
42
>>> translation = target_position - current_com
43
>>> protein.translate(translation)
44
"""
45
```
46
47
### Rotation
48
49
```python { .api }
50
def rotate(atomgroup, R, point=(0, 0, 0)):
51
"""
52
Rotate coordinates using a rotation matrix.
53
54
Parameters
55
----------
56
atomgroup : AtomGroup
57
Atoms to rotate.
58
R : array-like
59
3x3 rotation matrix. Must be orthogonal with determinant +1.
60
point : array-like, optional
61
Point about which to rotate (default origin).
62
63
Examples
64
--------
65
>>> import numpy as np
66
>>> # 90-degree rotation around z-axis
67
>>> R_z = np.array([[0, -1, 0],
68
... [1, 0, 0],
69
... [0, 0, 1]])
70
>>> protein.rotate(R_z)
71
72
>>> # Rotate around protein center
73
>>> center = protein.center_of_geometry()
74
>>> protein.rotate(R_z, point=center)
75
"""
76
77
def rotateby(atomgroup, angle, axis, point=None):
78
"""
79
Rotate by angle around axis.
80
81
Parameters
82
----------
83
atomgroup : AtomGroup
84
Atoms to rotate.
85
angle : float
86
Rotation angle in degrees.
87
axis : array-like
88
Rotation axis as 3-element unit vector.
89
point : array-like, optional
90
Point about which to rotate (default group center).
91
92
Examples
93
--------
94
>>> # Rotate 45 degrees around z-axis
95
>>> protein.rotateby(45.0, [0, 0, 1])
96
97
>>> # Rotate around custom axis
98
>>> axis = [1, 1, 0] / np.sqrt(2) # Normalize axis
99
>>> protein.rotateby(90.0, axis)
100
101
>>> # Rotate around specific point
102
>>> ca_atom = protein.select_atoms("name CA")[0]
103
>>> protein.rotateby(30.0, [0, 0, 1], point=ca_atom.position)
104
"""
105
```
106
107
### General Transformations
108
109
```python { .api }
110
def transform(atomgroup, M):
111
"""
112
Apply transformation matrix to coordinates.
113
114
Parameters
115
----------
116
atomgroup : AtomGroup
117
Atoms to transform.
118
M : array-like
119
Transformation matrix. Can be:
120
- 3x3 rotation/scaling matrix
121
- 4x4 homogeneous transformation matrix
122
- 3x4 [rotation|translation] matrix
123
124
Examples
125
--------
126
>>> # Combined rotation and translation using 4x4 matrix
127
>>> M = np.eye(4)
128
>>> M[:3, :3] = rotation_matrix # 3x3 rotation
129
>>> M[:3, 3] = translation_vector # Translation
130
>>> protein.transform(M)
131
132
>>> # Scaling transformation
133
>>> scale_matrix = np.diag([2.0, 1.0, 0.5]) # Scale x by 2, z by 0.5
134
>>> protein.transform(scale_matrix)
135
"""
136
```
137
138
## Structural Alignment
139
140
### Alignment Functions
141
142
```python { .api }
143
from MDAnalysis.analysis import align
144
145
def alignto(mobile, reference, select="all", weights=None, subsetA=None,
146
subsetB=None, **kwargs):
147
"""
148
Align mobile structure to reference structure.
149
150
Parameters
151
----------
152
mobile : Universe or AtomGroup
153
Structure to be aligned (modified in place).
154
reference : Universe or AtomGroup
155
Reference structure for alignment.
156
select : str, optional
157
Selection string for atoms used in alignment (default "all").
158
weights : str or array-like, optional
159
Weights for alignment. Can be "mass" for mass weighting or
160
array of weights matching selected atoms.
161
subsetA : AtomGroup, optional
162
Specific subset of mobile atoms for fitting.
163
subsetB : AtomGroup, optional
164
Specific subset of reference atoms for fitting.
165
**kwargs
166
Additional arguments for alignment algorithm.
167
168
Returns
169
-------
170
dict
171
Dictionary containing:
172
- 'rmsd': RMSD after alignment
173
- 'rotmat': 3x3 rotation matrix applied
174
- 'translation': translation vector applied
175
176
Examples
177
--------
178
>>> # Align protein backbones
179
>>> result = align.alignto(mobile_u, reference_u, select="backbone")
180
>>> print(f"Alignment RMSD: {result['rmsd']:.2f} Å")
181
182
>>> # Mass-weighted alignment of CA atoms
183
>>> result = align.alignto(mobile_u, reference_u,
184
... select="name CA", weights="mass")
185
186
>>> # Align using specific atom selections
187
>>> mobile_ca = mobile_u.select_atoms("name CA")
188
>>> ref_ca = reference_u.select_atoms("name CA")
189
>>> result = align.alignto(mobile_u, reference_u,
190
... subsetA=mobile_ca, subsetB=ref_ca)
191
"""
192
193
class AlignTraj:
194
"""
195
Align trajectory frames to reference structure.
196
197
Performs structural alignment of entire trajectory, optionally
198
writing aligned frames to output file.
199
"""
200
201
def __init__(self, mobile, reference, select="all", filename=None,
202
weights=None, in_memory=False, **kwargs):
203
"""
204
Initialize trajectory alignment.
205
206
Parameters
207
----------
208
mobile : Universe
209
Universe with trajectory to align.
210
reference : Universe or AtomGroup
211
Reference structure for alignment.
212
select : str, optional
213
Selection for alignment atoms (default "all").
214
filename : str, optional
215
Output filename for aligned trajectory.
216
weights : str or array-like, optional
217
Alignment weights ("mass" or custom array).
218
in_memory : bool, optional
219
Load trajectory in memory for speed (default False).
220
**kwargs
221
Additional alignment parameters.
222
223
Examples
224
--------
225
>>> # Align entire trajectory and save
226
>>> aligner = align.AlignTraj(u, reference,
227
... select="protein and name CA",
228
... filename="aligned.xtc")
229
>>> aligner.run()
230
231
>>> # In-memory alignment for analysis
232
>>> aligner = align.AlignTraj(u, reference,
233
... select="backbone",
234
... in_memory=True)
235
>>> aligner.run()
236
>>> rmsd_data = aligner.rmsd # RMSD per frame
237
"""
238
239
def run(self, start=None, stop=None, step=None, **kwargs):
240
"""
241
Execute trajectory alignment.
242
243
Parameters
244
----------
245
start : int, optional
246
First frame to align (default None for beginning).
247
stop : int, optional
248
Last frame to align (default None for end).
249
step : int, optional
250
Step size for alignment (default None for 1).
251
**kwargs
252
Additional run parameters.
253
254
Returns
255
-------
256
AlignTraj
257
Self for method chaining.
258
"""
259
260
@property
261
def rmsd(self):
262
"""
263
RMSD values for each aligned frame.
264
265
Returns
266
-------
267
numpy.ndarray
268
Array of RMSD values in Angstrom for each frame.
269
"""
270
```
271
272
### Rotation Matrix Calculations
273
274
```python { .api }
275
from MDAnalysis.analysis.align import rotation_matrix
276
277
def rotation_matrix(a, b, weights=None):
278
"""
279
Calculate optimal rotation matrix between two coordinate sets.
280
281
Uses Kabsch algorithm to find rotation matrix that minimizes
282
RMSD between coordinate sets.
283
284
Parameters
285
----------
286
a : array-like
287
Coordinate set A of shape (n, 3).
288
b : array-like
289
Coordinate set B of shape (n, 3).
290
weights : array-like, optional
291
Weights for each coordinate pair.
292
293
Returns
294
-------
295
tuple
296
(rotation_matrix, rmsd) where rotation_matrix is 3x3 array
297
and rmsd is root mean square deviation after alignment.
298
299
Examples
300
--------
301
>>> # Calculate rotation between CA atoms
302
>>> mobile_ca = mobile.select_atoms("name CA").positions
303
>>> ref_ca = reference.select_atoms("name CA").positions
304
>>> R, rmsd = rotation_matrix(mobile_ca, ref_ca)
305
>>> print(f"Optimal RMSD: {rmsd:.2f} Å")
306
307
>>> # Apply rotation to mobile structure
308
>>> mobile.atoms.rotate(R)
309
"""
310
311
def fit_rot_trans(mobile, reference, weights=None):
312
"""
313
Calculate optimal rotation and translation between structures.
314
315
Parameters
316
----------
317
mobile : array-like
318
Mobile coordinates of shape (n, 3).
319
reference : array-like
320
Reference coordinates of shape (n, 3).
321
weights : array-like, optional
322
Weights for fitting.
323
324
Returns
325
-------
326
tuple
327
(rotation, translation, rmsd) where rotation is 3x3 matrix,
328
translation is 3-element vector, and rmsd is final RMSD.
329
330
Examples
331
--------
332
>>> R, t, rmsd = fit_rot_trans(mobile_coords, ref_coords)
333
>>> # Apply transformation
334
>>> mobile_coords_aligned = mobile_coords @ R.T + t
335
"""
336
```
337
338
## Periodic Boundary Conditions
339
340
### Wrapping Coordinates
341
342
```python { .api }
343
def wrap(atomgroup, compound="atoms", center="com", box=None, inplace=True):
344
"""
345
Wrap coordinates into primary unit cell.
346
347
Parameters
348
----------
349
atomgroup : AtomGroup
350
Atoms to wrap.
351
compound : str, optional
352
Level for wrapping operation:
353
- "atoms": wrap individual atoms
354
- "group": wrap entire group as unit
355
- "residues": wrap by residue
356
- "segments": wrap by segment
357
- "molecules": wrap by molecule
358
- "fragments": wrap by fragment
359
center : str or array-like, optional
360
Center for wrapping:
361
- "com": center of mass
362
- "cog": center of geometry
363
- array: specific center point
364
box : array-like, optional
365
Unit cell dimensions. If None, uses current trajectory box.
366
inplace : bool, optional
367
Whether to modify coordinates in place (default True).
368
369
Returns
370
-------
371
numpy.ndarray or None
372
Wrapped coordinates if inplace=False, otherwise None.
373
374
Examples
375
--------
376
>>> # Wrap all atoms individually
377
>>> u.atoms.wrap(compound="atoms")
378
379
>>> # Wrap by residues to keep residues intact
380
>>> u.atoms.wrap(compound="residues", center="com")
381
382
>>> # Wrap around specific center
383
>>> center_point = [25.0, 25.0, 25.0]
384
>>> u.atoms.wrap(center=center_point)
385
386
>>> # Get wrapped coordinates without modifying original
387
>>> wrapped_coords = u.atoms.wrap(inplace=False)
388
"""
389
390
def unwrap(atomgroup, compound="fragments", reference="com", inplace=True):
391
"""
392
Unwrap coordinates across periodic boundaries.
393
394
Removes effects of periodic wrapping to make molecules whole
395
and trajectories continuous.
396
397
Parameters
398
----------
399
atomgroup : AtomGroup
400
Atoms to unwrap.
401
compound : str, optional
402
Level for unwrapping:
403
- "atoms": unwrap individual atoms
404
- "group": unwrap entire group
405
- "residues": unwrap by residue
406
- "segments": unwrap by segment
407
- "molecules": unwrap by molecule
408
- "fragments": unwrap by fragment (default)
409
reference : str or array-like, optional
410
Reference for unwrapping:
411
- "com": center of mass
412
- "cog": center of geometry
413
- array: specific reference point
414
inplace : bool, optional
415
Whether to modify coordinates in place (default True).
416
417
Returns
418
-------
419
numpy.ndarray or None
420
Unwrapped coordinates if inplace=False, otherwise None.
421
422
Examples
423
--------
424
>>> # Unwrap molecular fragments
425
>>> u.atoms.unwrap(compound="fragments")
426
427
>>> # Unwrap by residues
428
>>> protein = u.select_atoms("protein")
429
>>> protein.unwrap(compound="residues", reference="com")
430
431
>>> # Unwrap entire trajectory
432
>>> for ts in u.trajectory:
433
... u.atoms.unwrap(compound="molecules")
434
"""
435
436
def pack_into_box(atomgroup, box=None, inplace=True):
437
"""
438
Pack coordinates into simulation box.
439
440
Ensures all coordinates are within box boundaries by
441
applying appropriate translations.
442
443
Parameters
444
----------
445
atomgroup : AtomGroup
446
Atoms to pack into box.
447
box : array-like, optional
448
Box dimensions as [a, b, c, alpha, beta, gamma].
449
If None, uses current trajectory box.
450
inplace : bool, optional
451
Whether to modify coordinates in place (default True).
452
453
Returns
454
-------
455
numpy.ndarray or None
456
Packed coordinates if inplace=False, otherwise None.
457
458
Examples
459
--------
460
>>> # Pack all atoms into current box
461
>>> u.atoms.pack_into_box()
462
463
>>> # Pack into custom box
464
>>> box_dims = [50.0, 50.0, 50.0, 90.0, 90.0, 90.0]
465
>>> u.atoms.pack_into_box(box=box_dims)
466
"""
467
```
468
469
## On-the-Fly Transformations
470
471
### Transformation Classes
472
473
```python { .api }
474
from MDAnalysis import transformations
475
476
# Available transformation classes
477
class TransformationBase:
478
"""
479
Base class for trajectory transformations.
480
481
Transformations can be applied during trajectory reading
482
to modify coordinates on-the-fly.
483
"""
484
485
def __call__(self, ts):
486
"""
487
Apply transformation to timestep.
488
489
Parameters
490
----------
491
ts : Timestep
492
Timestep object to transform.
493
494
Returns
495
-------
496
Timestep
497
Modified timestep.
498
"""
499
500
class unwrap(TransformationBase):
501
"""Unwrap coordinates transformation."""
502
503
def __init__(self, atomgroup, **kwargs):
504
"""
505
Initialize unwrapping transformation.
506
507
Parameters
508
----------
509
atomgroup : AtomGroup
510
Atoms to unwrap in each frame.
511
**kwargs
512
Additional unwrapping parameters.
513
"""
514
515
class wrap(TransformationBase):
516
"""Wrap coordinates transformation."""
517
518
def __init__(self, atomgroup, **kwargs):
519
"""
520
Initialize wrapping transformation.
521
522
Parameters
523
----------
524
atomgroup : AtomGroup
525
Atoms to wrap in each frame.
526
**kwargs
527
Additional wrapping parameters.
528
"""
529
530
class center_in_box(TransformationBase):
531
"""Center system in simulation box."""
532
533
def __init__(self, atomgroup, center='mass', **kwargs):
534
"""
535
Initialize centering transformation.
536
537
Parameters
538
----------
539
atomgroup : AtomGroup
540
Atoms to center.
541
center : str, optional
542
Centering method ('mass' or 'geometry').
543
"""
544
545
class fit_rot_trans(TransformationBase):
546
"""Align trajectory to reference structure."""
547
548
def __init__(self, atomgroup, reference, **kwargs):
549
"""
550
Initialize alignment transformation.
551
552
Parameters
553
----------
554
atomgroup : AtomGroup
555
Atoms to align in each frame.
556
reference : AtomGroup or array-like
557
Reference coordinates for alignment.
558
"""
559
560
# Usage with Universe
561
transformations_list = [
562
transformations.unwrap(u.atoms),
563
transformations.center_in_box(u.select_atoms("protein")),
564
transformations.wrap(u.atoms, compound="residues")
565
]
566
567
u = mda.Universe("topology.psf", "trajectory.dcd",
568
transformations=transformations_list)
569
570
# Now trajectory is automatically transformed during reading
571
for ts in u.trajectory:
572
# Coordinates are already unwrapped, centered, and wrapped
573
pass
574
```
575
576
## Advanced Transformation Techniques
577
578
### Principal Axis Alignment
579
580
```python { .api }
581
def align_principal_axis(atomgroup, axis='z', vector=None):
582
"""
583
Align principal axis of molecule with coordinate axis.
584
585
Parameters
586
----------
587
atomgroup : AtomGroup
588
Atoms to align.
589
axis : str or int, optional
590
Target coordinate axis ('x', 'y', 'z' or 0, 1, 2).
591
vector : array-like, optional
592
Specific vector to align with (overrides axis).
593
594
Examples
595
--------
596
>>> # Align protein's principal axis with z-axis
597
>>> protein = u.select_atoms("protein")
598
>>> align_principal_axis(protein, axis='z')
599
600
>>> # Align with custom vector
601
>>> target_vector = [1, 1, 0] / np.sqrt(2)
602
>>> align_principal_axis(protein, vector=target_vector)
603
"""
604
# Calculate principal axes
605
principal_axes = atomgroup.principal_axes()
606
longest_axis = principal_axes[0] # Largest eigenvalue
607
608
if vector is None:
609
# Convert axis specification to vector
610
if axis == 'x' or axis == 0:
611
target = np.array([1, 0, 0])
612
elif axis == 'y' or axis == 1:
613
target = np.array([0, 1, 0])
614
elif axis == 'z' or axis == 2:
615
target = np.array([0, 0, 1])
616
else:
617
raise ValueError("axis must be 'x', 'y', 'z' or 0, 1, 2")
618
else:
619
target = np.array(vector)
620
target = target / np.linalg.norm(target)
621
622
# Calculate rotation matrix
623
from MDAnalysis.lib.mdamath import normal
624
rotation_axis = np.cross(longest_axis, target)
625
626
if np.allclose(rotation_axis, 0):
627
# Axes are already aligned or antiparallel
628
if np.dot(longest_axis, target) < 0:
629
# Antiparallel - rotate 180 degrees
630
perpendicular = normal(longest_axis)
631
angle = np.pi
632
else:
633
# Already aligned
634
return
635
else:
636
rotation_axis = rotation_axis / np.linalg.norm(rotation_axis)
637
angle = np.arccos(np.clip(np.dot(longest_axis, target), -1, 1))
638
639
# Create rotation matrix using Rodrigues' formula
640
cos_angle = np.cos(angle)
641
sin_angle = np.sin(angle)
642
643
# Cross product matrix
644
K = np.array([[0, -rotation_axis[2], rotation_axis[1]],
645
[rotation_axis[2], 0, -rotation_axis[0]],
646
[-rotation_axis[1], rotation_axis[0], 0]])
647
648
R = np.eye(3) + sin_angle * K + (1 - cos_angle) * np.dot(K, K)
649
650
# Apply rotation
651
center = atomgroup.center_of_geometry()
652
atomgroup.translate(-center)
653
atomgroup.rotate(R)
654
atomgroup.translate(center)
655
```
656
657
### Coordinate System Transformations
658
659
```python { .api }
660
def cartesian_to_spherical(coordinates, origin=None):
661
"""
662
Convert Cartesian to spherical coordinates.
663
664
Parameters
665
----------
666
coordinates : array-like
667
Cartesian coordinates of shape (..., 3).
668
origin : array-like, optional
669
Origin for coordinate system (default [0, 0, 0]).
670
671
Returns
672
-------
673
numpy.ndarray
674
Spherical coordinates [r, theta, phi] where:
675
- r: radial distance
676
- theta: polar angle (0 to π)
677
- phi: azimuthal angle (0 to 2π)
678
679
Examples
680
--------
681
>>> coords = protein.positions
682
>>> com = protein.center_of_mass()
683
>>> spherical = cartesian_to_spherical(coords, origin=com)
684
>>> r = spherical[..., 0] # Distances from center
685
>>> theta = spherical[..., 1] # Polar angles
686
>>> phi = spherical[..., 2] # Azimuthal angles
687
"""
688
coords = np.asarray(coordinates)
689
if origin is not None:
690
coords = coords - np.asarray(origin)
691
692
# Radial distance
693
r = np.sqrt(np.sum(coords**2, axis=-1))
694
695
# Polar angle (theta)
696
theta = np.arccos(coords[..., 2] / r)
697
698
# Azimuthal angle (phi)
699
phi = np.arctan2(coords[..., 1], coords[..., 0])
700
phi = np.where(phi < 0, phi + 2*np.pi, phi) # Ensure 0 to 2π
701
702
return np.stack([r, theta, phi], axis=-1)
703
704
def cylindrical_coordinates(coordinates, axis='z', origin=None):
705
"""
706
Convert to cylindrical coordinates.
707
708
Parameters
709
----------
710
coordinates : array-like
711
Cartesian coordinates of shape (..., 3).
712
axis : str or int, optional
713
Cylinder axis ('x', 'y', 'z' or 0, 1, 2).
714
origin : array-like, optional
715
Origin for coordinate system.
716
717
Returns
718
-------
719
numpy.ndarray
720
Cylindrical coordinates [rho, phi, z] where:
721
- rho: radial distance from axis
722
- phi: azimuthal angle
723
- z: position along axis
724
725
Examples
726
--------
727
>>> # Analyze protein in cylindrical coordinates around z-axis
728
>>> coords = protein.positions
729
>>> cylindrical = cylindrical_coordinates(coords, axis='z')
730
>>> rho = cylindrical[..., 0] # Distance from z-axis
731
>>> phi = cylindrical[..., 1] # Angle around z-axis
732
>>> z = cylindrical[..., 2] # Height along z-axis
733
"""
734
coords = np.asarray(coordinates)
735
if origin is not None:
736
coords = coords - np.asarray(origin)
737
738
# Determine axis index
739
if axis == 'x' or axis == 0:
740
ax_idx = 0
741
perp_indices = [1, 2]
742
elif axis == 'y' or axis == 1:
743
ax_idx = 1
744
perp_indices = [0, 2]
745
elif axis == 'z' or axis == 2:
746
ax_idx = 2
747
perp_indices = [0, 1]
748
else:
749
raise ValueError("axis must be 'x', 'y', 'z' or 0, 1, 2")
750
751
# Axial coordinate
752
z = coords[..., ax_idx]
753
754
# Radial distance in perpendicular plane
755
rho = np.sqrt(coords[..., perp_indices[0]]**2 +
756
coords[..., perp_indices[1]]**2)
757
758
# Azimuthal angle
759
phi = np.arctan2(coords[..., perp_indices[1]],
760
coords[..., perp_indices[0]])
761
phi = np.where(phi < 0, phi + 2*np.pi, phi)
762
763
return np.stack([rho, phi, z], axis=-1)
764
```
765
766
## Usage Examples
767
768
### Complete Alignment Workflow
769
770
```python { .api }
771
def alignment_workflow(mobile_universe, reference_universe):
772
"""
773
Complete structural alignment workflow example.
774
"""
775
# Step 1: Initial alignment using backbone atoms
776
print("Performing initial backbone alignment...")
777
result = align.alignto(mobile_universe, reference_universe,
778
select="backbone")
779
print(f"Initial RMSD: {result['rmsd']:.2f} Å")
780
781
# Step 2: Refined alignment using CA atoms with mass weighting
782
print("Refining alignment with CA atoms...")
783
result = align.alignto(mobile_universe, reference_universe,
784
select="name CA", weights="mass")
785
print(f"Refined RMSD: {result['rmsd']:.2f} Å")
786
787
# Step 3: Calculate RMSD for different regions
788
regions = {
789
'backbone': 'backbone',
790
'sidechains': 'protein and not backbone',
791
'alpha_carbons': 'name CA',
792
'binding_site': 'resid 23 45 67 89'
793
}
794
795
print("\nRegion-specific RMSDs after alignment:")
796
for region_name, selection in regions.items():
797
mobile_sel = mobile_universe.select_atoms(selection)
798
ref_sel = reference_universe.select_atoms(selection)
799
800
if len(mobile_sel) > 0 and len(ref_sel) > 0:
801
# Calculate RMSD manually
802
diff = mobile_sel.positions - ref_sel.positions
803
rmsd = np.sqrt(np.mean(np.sum(diff**2, axis=1)))
804
print(f" {region_name}: {rmsd:.2f} Å")
805
806
return result
807
808
# Run alignment workflow
809
result = alignment_workflow(mobile_u, reference_u)
810
```
811
812
### Trajectory Preprocessing Pipeline
813
814
```python { .api }
815
def preprocess_trajectory(universe, output_file):
816
"""
817
Comprehensive trajectory preprocessing pipeline.
818
"""
819
protein = universe.select_atoms("protein")
820
821
# Define transformation pipeline
822
transformations_list = [
823
# 1. Unwrap molecules to make them whole
824
transformations.unwrap(universe.atoms, compound="molecules"),
825
826
# 2. Center protein in box
827
transformations.center_in_box(protein, center='mass'),
828
829
# 3. Align protein backbone to reference (first frame)
830
transformations.fit_rot_trans(
831
protein.select_atoms("backbone"),
832
protein.select_atoms("backbone").positions # First frame as reference
833
),
834
835
# 4. Wrap solvent around protein
836
transformations.wrap(universe.atoms, compound="residues")
837
]
838
839
# Apply transformations and write processed trajectory
840
with mda.Writer(output_file, n_atoms=universe.atoms.n_atoms) as W:
841
for ts in universe.trajectory:
842
# Apply transformations
843
for transform in transformations_list:
844
ts = transform(ts)
845
846
# Write processed frame
847
W.write(universe.atoms, ts=ts)
848
849
print(f"Processed trajectory written to: {output_file}")
850
851
# Process trajectory
852
preprocess_trajectory(u, "processed_trajectory.xtc")
853
```
854
855
### Periodic Boundary Analysis
856
857
```python { .api }
858
def analyze_pbc_effects(universe):
859
"""
860
Analyze effects of periodic boundary conditions.
861
"""
862
protein = universe.select_atoms("protein")
863
864
results = {
865
'frame': [],
866
'original_rgyr': [],
867
'wrapped_rgyr': [],
868
'unwrapped_rgyr': [],
869
'box_violations': []
870
}
871
872
for ts in universe.trajectory:
873
frame = ts.frame
874
box_dims = ts.dimensions[:3] # Box lengths
875
876
# Original coordinates
877
original_rgyr = protein.radius_of_gyration()
878
879
# Wrapped coordinates
880
protein_wrapped = protein.copy()
881
protein_wrapped.wrap(compound="residues")
882
wrapped_rgyr = protein_wrapped.radius_of_gyration()
883
884
# Unwrapped coordinates
885
protein_unwrapped = protein.copy()
886
protein_unwrapped.unwrap(compound="residues")
887
unwrapped_rgyr = protein_unwrapped.radius_of_gyration()
888
889
# Check for atoms outside box
890
coords = protein.positions
891
violations = 0
892
for dim in range(3):
893
if box_dims[dim] > 0: # Valid box dimension
894
outside = np.sum((coords[:, dim] < 0) |
895
(coords[:, dim] > box_dims[dim]))
896
violations += outside
897
898
# Store results
899
results['frame'].append(frame)
900
results['original_rgyr'].append(original_rgyr)
901
results['wrapped_rgyr'].append(wrapped_rgyr)
902
results['unwrapped_rgyr'].append(unwrapped_rgyr)
903
results['box_violations'].append(violations)
904
905
# Analysis summary
906
print("Periodic Boundary Condition Analysis:")
907
print(f"Average Rgyr (original): {np.mean(results['original_rgyr']):.2f} Å")
908
print(f"Average Rgyr (wrapped): {np.mean(results['wrapped_rgyr']):.2f} Å")
909
print(f"Average Rgyr (unwrapped): {np.mean(results['unwrapped_rgyr']):.2f} Å")
910
print(f"Frames with box violations: {np.sum(np.array(results['box_violations']) > 0)}")
911
912
return results
913
914
# Analyze PBC effects
915
pbc_analysis = analyze_pbc_effects(u)
916
```
917
918
### Custom Transformation Development
919
920
```python { .api }
921
class RemoveCOMMotion(transformations.TransformationBase):
922
"""
923
Custom transformation to remove center-of-mass motion.
924
925
Removes translational motion by centering system at origin
926
and removes rotational motion by aligning to reference.
927
"""
928
929
def __init__(self, atomgroup, reference_frame=0):
930
"""
931
Initialize COM motion removal.
932
933
Parameters
934
----------
935
atomgroup : AtomGroup
936
Atoms for COM calculation and motion removal.
937
reference_frame : int, optional
938
Frame to use as reference for rotation removal.
939
"""
940
self.atomgroup = atomgroup
941
self.reference_frame = reference_frame
942
self.reference_coords = None
943
self._first_frame = True
944
945
def __call__(self, ts):
946
"""Apply transformation to timestep."""
947
# Get current frame coordinates
948
current_coords = self.atomgroup.positions.copy()
949
950
# Remove translational motion (center at origin)
951
com = np.mean(current_coords, axis=0)
952
current_coords -= com
953
954
# Set up reference for rotational alignment
955
if self._first_frame and ts.frame == self.reference_frame:
956
self.reference_coords = current_coords.copy()
957
self._first_frame = False
958
959
# Remove rotational motion (align to reference)
960
if self.reference_coords is not None:
961
from MDAnalysis.analysis.align import rotation_matrix
962
R, _ = rotation_matrix(current_coords, self.reference_coords)
963
current_coords = current_coords @ R.T
964
965
# Update positions in timestep
966
self.atomgroup.positions = current_coords
967
968
return ts
969
970
# Use custom transformation
971
protein = u.select_atoms("protein and name CA")
972
custom_transform = RemoveCOMMotion(protein, reference_frame=0)
973
974
# Apply during trajectory reading
975
u_processed = mda.Universe("topology.psf", "trajectory.dcd",
976
transformations=[custom_transform])
977
978
for ts in u_processed.trajectory:
979
# Protein motion has been removed
980
pass
981
```