Examples
Ready-to-run example scripts demonstrating common workflows.
Parsing Trajectory Files
This example demonstrates how to parse different trajectory file formats.
1"""
2Example: Using DumpParser and DumpWaterMoleculeFinder
3
4This example shows how to:
51. Identify water molecules in a LAMMPS dump file.
62. Extract only their oxygen atom coordinates.
7"""
8
9from wetting_angle_kit.parser import (
10 ASE_Parser,
11 ASE_WaterMoleculeFinder,
12 DumpParser,
13 DumpWaterMoleculeFinder,
14 XYZParser,
15)
16
17# --- Define input file ---
18filename = "../../tests/trajectories/traj_10_3_330w_nve_4k_reajust.lammpstrj"
19
20# --- Initialize water molecule finder ---
21wat_find = DumpWaterMoleculeFinder(
22 filename,
23 particle_type_wall={3}, # atom type for wall
24 oxygen_type=1, # atom type for oxygen
25 hydrogen_type=2, # atom type for hydrogen
26)
27
28# --- Identify water oxygen indices for the first frame ---
29oxygen_indices = wat_find.get_water_oxygen_ids(frame_index=0)
30print(f"Number of water molecules: {len(oxygen_indices)}")
31
32# --- Initialize parser ---
33parser = DumpParser(filename)
34
35# --- Extract only oxygen coordinates for frame 0 ---
36oxygen_positions = parser.parse(frame_index=0, indices=oxygen_indices)
37print("Extracted oxygen coordinates shape:", oxygen_positions.shape)
38
39# --- Optional: Extract all atoms ---
40# all_positions = parser.parse(frame_index=0)
41# print("All atom positions shape:", all_positions.shape)
42
43"""
44Example: Using ASE_Parser and ASE_WaterMoleculeFinder
45
46This example demonstrates how to:
471. Identify water oxygens in an ASE trajectory.
482. Extract their positions for a given frame.
49"""
50
51# --- Define input file ---
52filename = "../../tests/trajectories/slice_10_mace_mlips_cylindrical_2_5.traj"
53
54# --- Initialize water molecule finder ---
55wat_find = ASE_WaterMoleculeFinder(
56 filename,
57 particle_type_wall=["C"], # element name for wall
58 oh_cutoff=0.4, # O–H cutoff distance (Å)
59)
60
61# --- Get oxygen indices for frame 0 ---
62oxygen_indices = wat_find.get_water_oxygen_indices(frame_index=0)
63print(f"Number of water molecules: {len(oxygen_indices)}")
64
65# --- Initialize parser ---
66parser = ASE_Parser(filename)
67
68# --- Extract oxygen coordinates only ---
69oxygen_positions = parser.parse(frame_index=0, indices=oxygen_indices)
70print("Extracted oxygen coordinates shape:", oxygen_positions.shape)
71
72"""
73Example: Using XYZParser
74
75This example demonstrates how to:
761. Load atomic positions from an XYZ file.
772. Extract all atoms or a subset of atoms.
78"""
79
80# --- Define input file ---
81filename = "../../tests/trajectories/slice_10_mace_mlips_cylindrical_2_5.xyz"
82
83# --- Initialize parser ---
84xyz_parser = XYZParser(filename)
85
86# --- Extract all atom coordinates for frame 0 ---
87positions = xyz_parser.parse(frame_index=0)
88print("Total atoms loaded:", len(positions))
89
90# --- Extract subset of atoms (first 50) ---
91subset = xyz_parser.parse(frame_index=0, indices=list(range(50)))
92print("Subset (50 atoms) shape:", subset.shape)
Binning Contact Angle Analysis
Example of using the binning method for contact angle analysis.
1# Import necessary modules
2from wetting_angle_kit.contact_angle_method import contact_angle_analyzer
3from wetting_angle_kit.parser import DumpParser, DumpWaterMoleculeFinder
4
5# --- Step 1: Define the trajectory file ---
6filename = "../../tests/trajectories/traj_spherical_drop_4k.lammpstrj"
7
8# --- Step 2: Initialize the water molecule finder ---
9# This identifies O and H atoms in water molecules
10wat_find = DumpWaterMoleculeFinder(
11 filename,
12 particle_type_wall={3}, # Wall atom types
13 oxygen_type=1, # Oxygen atom type
14 hydrogen_type=2, # Hydrogen atom type
15)
16
17# --- Step 3: Get oxygen atom indices for the first frame ---
18oxygen_indices = wat_find.get_water_oxygen_ids(frame_index=0)
19print("Number of water molecules:", len(oxygen_indices))
20
21# --- Step 4: Define binning parameters ---
22binning_params = {
23 "xi_0": 0.0, # Minimum x-coordinate
24 "xi_f": 70.0, # Maximum x-coordinate
25 "nbins_xi": 30, # Number of bins along x
26 "zi_0": 0.0, # Minimum z-coordinate
27 "zi_f": 70.0, # Maximum z-coordinate
28 "nbins_zi": 30, # Number of bins along z
29}
30
31# --- Step 5: Initialize the parser ---
32parser = DumpParser(filename)
33
34# --- Step 6: Create the contact angle analyzer ---
35analyzer = contact_angle_analyzer(
36 method="binning",
37 parser=parser,
38 output_dir="results_binning_example",
39 atom_indices=oxygen_indices,
40 droplet_geometry="spherical", # Interface fitting model
41 binning_params=binning_params,
42 plot_graphs=True, # Enable plotting for automated runs
43)
44
45# --- Step 7: Run analysis for a frame range ---
46results = analyzer.analyze([1]) # Analyze frame 1
47print("Analysis results:", results)
Sliced Contact Angle Analysis
Example of using the sliced method for contact angle analysis.
1
Visualizing Sliced Trajectories
Example of visualizing droplet trajectories with the sliced method.
1import numpy as np
2
3from wetting_angle_kit.contact_angle_method.sliced_method import ContactAngleSliced
4from wetting_angle_kit.parser import (
5 DumpParser,
6 DumpWallParser,
7 DumpWaterMoleculeFinder,
8)
9from wetting_angle_kit.visualization_statistics_angles import Droplet_sliced_Plotter
10
11# --- 1. Define the Input Trajectory ---
12# Note: Ensure this path points to your actual .lammpstrj file location
13filename = (
14 "../wetting_angle_kit/tests/trajectories/traj_10_3_330w_nve_4k_reajust.lammpstrj"
15)
16
17# --- 2. Identify Water Molecules ---
18wat_find = DumpWaterMoleculeFinder(
19 filename, particle_type_wall={3}, oxygen_type=1, hydrogen_type=2
20)
21
22oxygen_indices = wat_find.get_water_oxygen_ids(frame_index=0)
23print("Number of water molecules detected:", len(oxygen_indices))
24
25# --- 3. Parse Atomic Coordinates ---
26parser = DumpParser(filepath=filename)
27oxygen_position = parser.parse(frame_index=10, indices=oxygen_indices)
28
29coord_wall = DumpWallParser(filename, particule_liquid_type={1, 2})
30wall_coords = coord_wall.parse(frame_index=1)
31
32# --- 4. Compute Contact Angles ---
33#
34
35
36processor = ContactAngleSliced(
37 liquid_coordinates=oxygen_position,
38 liquid_geom_center=np.mean(oxygen_position, axis=0),
39 droplet_geometry="cylinder_y",
40 delta_cylinder=5,
41 max_dist=100,
42 width_cylinder=21,
43)
44
45list_alfas, array_surfaces, array_popt = processor.predict_contact_angle()
46print("Mean contact angles (°):", list_alfas)
47
48# --- 5. Visualize the Droplet ---
49plotter = Droplet_sliced_Plotter(center=True, show_wall=True, molecule_view=True)
50
51plotter.plot_surface_points(
52 oxygen_position=oxygen_position,
53 surface_data=array_surfaces,
54 popt=array_popt[0],
55 wall_coords=wall_coords,
56 output_filename="droplet_plot.png",
57 alpha=list_alfas[0],
58)
59
60print(" Plot saved as 'droplet_plot.png'")