RDF (Radial Distribution Function)#

Radial Distribution Functions are crucial to represent the structural correlations in the AIMD trajectories accurately. We reconstructed the results of the original work by Villard et al. [1] by computing the RDFs for all three pair types (O–O, O–H, H–H) across the five meta-GGA functionals. Our CRISP results are presented alongside experimental references at 298K provided by the original work gathered via X-ray diffraction (for g_OO), interpolated at 298K and joint X-ray/neutron diffraction experiments (for g_OH and g_HH).

Basic RDF Analysis#

from CRISP.data_analysis.prdf import analyze_rdf
import numpy as np
import os

# Minnesota functionals from Villard et al. study
functionals = ['M06-L', 'M11-L', 'MN12-L', 'MN15-L', 'revM06-L']

# Directory mapping for file paths
func_dirs = {
    'M06-L': 'M06L',
    'M11-L': 'M11L',
    'MN12-L': 'MN12L',
    'MN15-L': 'MN15L',
    'revM06-L': 'REVM06L'
}

peak_data = {}

print("Analyzing O-H radial distribution functions...")

for functional in functionals:
    dir_name = func_dirs[functional]
    print(f"\nWorking on {functional}...")

    # Set up file paths
    oxygen_indices = f'./Data_supplementary_analysis/MetaGGA/{dir_name}/indices_new/O_indices.npy'
    hydrogen_indices = f'./Data_supplementary_analysis/MetaGGA/{dir_name}/indices_new/H_indices.npy'
    trajectory = f'./supplementary_data/MetaGGA/{dir_name}/TRAJEC.traj'
    results_dir = f'./Data_supplementary_analysis/MetaGGA/{dir_name}/goh'

    os.makedirs(results_dir, exist_ok=True)

    try:
        # Load atom indices
        o_atoms = np.load(oxygen_indices)
        h_atoms = np.load(hydrogen_indices)

        # Setup for O-H partial RDF calculation
        atom_pairs = (o_atoms.tolist(), h_atoms.tolist())

        # Calculate RDF
        rdf_results = analyze_rdf(
            use_prdf=True,
            rmax=6.0,
            traj_path=trajectory,
            nbins=50,
            frame_skip=1,
            output_filename="prdf_o-atoms_h-atoms",
            atomic_indices=atom_pairs,
            output_dir=results_dir,
            create_plots=True
        )

        # Extract first coordination shell peak
        distances = rdf_results['x_data']
        avg_rdf = np.mean(rdf_results['y_data_all'], axis=0)
        first_peak_idx = np.argmax(avg_rdf)
        first_peak_pos = distances[first_peak_idx]

        peak_data[functional] = first_peak_pos
        print(f"  First O-H peak: {first_peak_pos:.2f} Å")

    except Exception as error:
        print(f"  Failed to process {functional}: {error}")

# Results summary
print("\n" + "="*40)
print("O-H COORDINATION ANALYSIS SUMMARY")
print("="*40)
print("Functional\t\tFirst Peak (Å)")
print("-"*40)

for func in functionals:
    if func in peak_data:
        tab = "\t\t" if len(func) <= 6 else "\t"
        print(f"{func}{tab}{peak_data[func]:.2f}")

print(f"\nAnalyzed {len(peak_data)} functionals successfully")

Similarly, we calculated gOO and gHH by modifying the atom indices appropriately:

  • gOO analysis: Used oxygen indices for both reference and target atoms

  • gHH analysis: Used hydrogen indices for both reference and target atoms

Results and Comparison#

O–O RDFs: CRISP replicates the overall features reported by Villard et al. [1], with clear functional differences. M06-L and M11-L understructure water (weaker first peak, higher first minimum), while MN12-L and MN15-L lead to stronger, slightly overstructured networks. revM06-L achieves intermediate structuring. CRISP-derived RDFs systematically shift O–O peaks slightly closer to experimental positions.

O–H and H–H RDFs: CRISP and reported values are almost similar for O-H and H-H peaks, with minor details such as revM06-L showing more pronounced second peak (around 4.0 Å) aligning with experimental trends.

Visualization:

Radial distribution functions computed using CRISP

Summary#

CRISP demonstrates optimized binning protocols, automated handling of periodic boundary conditions and statistical averaging, providing accurate reproduction of experimental trends across different theoretical methods.

References#