omecp/
pes_scan.rs

1//! PES (Potential Energy Surface) scanning functionality.
2//!
3//! This module implements comprehensive PES scanning capabilities following
4//!
5//! - 1D and 2D potential energy surface scans
6//! - Automatic scan grid generation
7//! - Constrained MECP optimization at each scan point
8//! - Result analysis and visualization data generation
9//! - Comprehensive reporting and convergence tracking
10
11use crate::geometry::Geometry;
12
13/// Represents the result of a single PES scan point.
14#[derive(Debug, Clone)]
15pub struct ScanPointResult {
16    /// First scan coordinate value
17    pub coord1: f64,
18    /// Second scan coordinate value (may be dummy for 1D scans)
19    pub coord2: f64,
20    /// Energy of state A at this point
21    pub energy_a: f64,
22    /// Energy of state B at this point
23    pub energy_b: f64,
24    /// Energy difference (E_A - E_B)
25    pub energy_diff: f64,
26    /// Whether the optimization converged at this point
27    pub converged: bool,
28    /// Number of optimization steps taken
29    pub num_steps: usize,
30    /// Final geometry at this scan point
31    pub geometry: Geometry,
32}
33
34/// Analyzes and collects PES scan results.
35///
36/// This function implements scan result analysis capabilities including:
37/// - Energy surface data collection
38/// - Convergence tracking for each scan point
39/// - Summary report generation
40/// - Data formatting for plotting
41///
42/// # Arguments
43///
44/// * `scan_results` - Collection of scan point results
45/// * `output_file` - Path for the analysis output file
46///
47/// # Returns
48///
49/// Returns `Ok(())` on successful analysis completion.
50pub fn analyze_scan_results(
51    scan_results: &[ScanPointResult],
52    output_file: &str,
53) -> Result<(), Box<dyn std::error::Error>> {
54    println!("\n****Analyzing PES Scan Results****");
55
56    // Generate scan summary report
57    generate_scan_summary(scan_results, output_file)?;
58
59    // Generate energy surface data for plotting
60    generate_energy_surface_data(scan_results, &format!("{}.dat", output_file))?;
61
62    // Generate convergence tracking report
63    generate_convergence_report(scan_results, &format!("{}_convergence.txt", output_file))?;
64
65    println!("Scan analysis completed. Results saved to: {}", output_file);
66    Ok(())
67}
68
69/// Generates a comprehensive scan summary report.
70///
71/// This function creates a detailed summary of the PES scan results,
72/// including statistics, energy ranges, and convergence information.
73///
74/// # Arguments
75///
76/// * `scan_results` - Collection of scan point results
77/// * `output_file` - Path for the summary report
78///
79/// # Returns
80///
81/// Returns `Ok(())` on successful report generation.
82fn generate_scan_summary(
83    scan_results: &[ScanPointResult],
84    output_file: &str,
85) -> Result<(), Box<dyn std::error::Error>> {
86    let mut summary = String::new();
87
88    summary.push_str("# PES Scan Summary Report\n");
89    summary.push_str("# Generated by OpenMECP\n\n");
90
91    // Basic statistics
92    summary.push_str(&format!("Total scan points: {}\n", scan_results.len()));
93
94    let converged_count = scan_results.iter().filter(|r| r.converged).count();
95    summary.push_str(&format!("Converged points: {}\n", converged_count));
96    summary.push_str(&format!(
97        "Convergence rate: {:.1}%\n\n",
98        100.0 * converged_count as f64 / scan_results.len() as f64
99    ));
100
101    // Energy statistics (only for converged points)
102    let converged_results: Vec<_> = scan_results.iter().filter(|r| r.converged).collect();
103    if !converged_results.is_empty() {
104        let energies_a: Vec<f64> = converged_results.iter().map(|r| r.energy_a).collect();
105        let energies_b: Vec<f64> = converged_results.iter().map(|r| r.energy_b).collect();
106        let energy_diffs: Vec<f64> = converged_results.iter().map(|r| r.energy_diff).collect();
107
108        let min_a = energies_a.iter().fold(f64::INFINITY, |a, &b| a.min(b));
109        let max_a = energies_a.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
110        let min_b = energies_b.iter().fold(f64::INFINITY, |a, &b| a.min(b));
111        let max_b = energies_b.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
112        let min_diff = energy_diffs.iter().fold(f64::INFINITY, |a, &b| a.min(b));
113        let max_diff = energy_diffs
114            .iter()
115            .fold(f64::NEG_INFINITY, |a, &b| a.max(b));
116
117        summary.push_str("Energy Statistics (Hartree, converged points only):\n");
118        summary.push_str(&format!("State A: {:.6} to {:.6}\n", min_a, max_a));
119        summary.push_str(&format!("State B: {:.6} to {:.6}\n", min_b, max_b));
120        summary.push_str(&format!(
121            "Energy Difference: {:.6} to {:.6}\n\n",
122            min_diff, max_diff
123        ));
124
125        // Find minimum energy difference point
126        if let Some(min_diff_result) = converged_results.iter().min_by(|a, b| {
127            a.energy_diff
128                .abs()
129                .partial_cmp(&b.energy_diff.abs())
130                .unwrap()
131        }) {
132            summary.push_str("Minimum |ΔE| Point:\n");
133            summary.push_str(&format!(
134                "Coordinates: ({:.4}, {:.4})\n",
135                min_diff_result.coord1, min_diff_result.coord2
136            ));
137            summary.push_str(&format!(
138                "Energy difference: {:.6} Hartree ({:.3} eV)\n\n",
139                min_diff_result.energy_diff,
140                min_diff_result.energy_diff * 27.211386
141            ));
142        }
143    }
144
145    // Detailed point-by-point results
146    summary.push_str("Detailed Results:\n");
147    summary.push_str("Coord1    Coord2    Energy_A      Energy_B      Delta_E       Conv  Steps\n");
148    summary.push_str("------------------------------------------------------------------------\n");
149
150    for result in scan_results {
151        if result.converged {
152            summary.push_str(&format!(
153                "{:8.4} {:8.4} {:12.6} {:12.6} {:12.6} {:4} {:5}\n",
154                result.coord1,
155                result.coord2,
156                result.energy_a,
157                result.energy_b,
158                result.energy_diff,
159                "Yes",
160                result.num_steps
161            ));
162        } else {
163            summary.push_str(&format!(
164                "{:8.4} {:8.4} {:>12} {:>12} {:>12} {:4} {:5}\n",
165                result.coord1, result.coord2, "FAILED", "FAILED", "FAILED", "No", result.num_steps
166            ));
167        }
168    }
169
170    std::fs::write(output_file, summary)?;
171    println!("Scan summary saved to: {}", output_file);
172    Ok(())
173}
174
175/// Generates energy surface data suitable for plotting.
176///
177/// This function creates data files in formats suitable for plotting programs
178/// like gnuplot, matplotlib, or other visualization tools.
179///
180/// # Arguments
181///
182/// * `scan_results` - Collection of scan point results
183/// * `output_file` - Path for the plotting data file
184///
185/// # Returns
186///
187/// Returns `Ok(())` on successful data file generation.
188fn generate_energy_surface_data(
189    scan_results: &[ScanPointResult],
190    output_file: &str,
191) -> Result<(), Box<dyn std::error::Error>> {
192    let mut data = String::new();
193
194    // Header for plotting data
195    data.push_str("# PES Scan Energy Surface Data\n");
196    data.push_str("# Columns: Coord1 Coord2 Energy_A Energy_B Energy_Diff Converged\n");
197    data.push_str("# Suitable for gnuplot splot or matplotlib\n");
198    data.push_str("# Non-converged points have energy values set to NaN\n\n");
199
200    for result in scan_results {
201        if result.converged {
202            data.push_str(&format!(
203                "{:.6} {:.6} {:.8} {:.8} {:.8} 1\n",
204                result.coord1, result.coord2, result.energy_a, result.energy_b, result.energy_diff
205            ));
206        } else {
207            data.push_str(&format!(
208                "{:.6} {:.6} NaN NaN NaN 0\n",
209                result.coord1, result.coord2
210            ));
211        }
212    }
213
214    std::fs::write(output_file, data)?;
215    println!("Energy surface data saved to: {}", output_file);
216    Ok(())
217}
218
219/// Generates a convergence tracking report.
220///
221/// This function analyzes the convergence behavior across the scan,
222/// identifying problematic regions and optimization statistics.
223///
224/// # Arguments
225///
226/// * `scan_results` - Collection of scan point results
227/// * `output_file` - Path for the convergence report
228///
229/// # Returns
230///
231/// Returns `Ok(())` on successful report generation.
232fn generate_convergence_report(
233    scan_results: &[ScanPointResult],
234    output_file: &str,
235) -> Result<(), Box<dyn std::error::Error>> {
236    let mut report = String::new();
237
238    report.push_str("# PES Scan Convergence Analysis\n\n");
239
240    // Convergence statistics
241    let total_points = scan_results.len();
242    let converged_points = scan_results.iter().filter(|r| r.converged).count();
243    let failed_points = total_points - converged_points;
244
245    report.push_str(&format!("Total scan points: {}\n", total_points));
246    report.push_str(&format!("Successfully converged: {}\n", converged_points));
247    report.push_str(&format!("Failed to converge: {}\n", failed_points));
248    report.push_str(&format!(
249        "Success rate: {:.1}%\n\n",
250        100.0 * converged_points as f64 / total_points as f64
251    ));
252
253    // Step count statistics (only for converged points)
254    let converged_results: Vec<_> = scan_results.iter().filter(|r| r.converged).collect();
255    if !converged_results.is_empty() {
256        let step_counts: Vec<usize> = converged_results.iter().map(|r| r.num_steps).collect();
257        let avg_steps = step_counts.iter().sum::<usize>() as f64 / step_counts.len() as f64;
258        let min_steps = *step_counts.iter().min().unwrap_or(&0);
259        let max_steps = *step_counts.iter().max().unwrap_or(&0);
260
261        report.push_str("Optimization Steps Statistics (converged points only):\n");
262        report.push_str(&format!("Average steps: {:.1}\n", avg_steps));
263        report.push_str(&format!("Minimum steps: {}\n", min_steps));
264        report.push_str(&format!("Maximum steps: {}\n\n", max_steps));
265    }
266
267    // List problematic points
268    report.push_str("Non-converged Points:\n");
269    if failed_points > 0 {
270        report.push_str("Coord1    Coord2    Steps\n");
271        report.push_str("------------------------\n");
272        for result in scan_results.iter().filter(|r| !r.converged) {
273            report.push_str(&format!(
274                "{:8.4} {:8.4} {:5}\n",
275                result.coord1, result.coord2, result.num_steps
276            ));
277        }
278    } else {
279        report.push_str("All scan points converged successfully!\n");
280    }
281
282    std::fs::write(output_file, report)?;
283    println!("Convergence report saved to: {}", output_file);
284    Ok(())
285}
omecp/pes_scan.rs

omecp/
pes_scan.rs
