README.md

Fortran 2D Wavepacket Propagation with Diabatic Surfaces and Berry Curvature

Overview

This program simulates the time-dependent quantum dynamics of a 2D wavepacket on two coupled diabatic potential energy surfaces. It employs the Crank-Nicolson propagation scheme (split-operator method with kinetic and potential propagators). The simulation calculates various observables, including:

• Population dynamics on diabatic states
• Transmission probabilities into different product channels
• Spin selectivity (based on channel populations and flux)
• Berry curvature of the lower adiabatic potential energy surface

The code is a Fortran conversion of an original Python script, designed for potentially improved performance for large-scale simulations. It utilizes OpenMP for parallelization over different initial collision energies. All plotting is deferred; instead, data is saved to CSV files for later analysis and visualization.

Code Structure

The project is organized into several Fortran modules, each responsible for specific functionalities:

Module	Description
main_crank_prop.f90	The main program that orchestrates the simulation. It handles the overall setup, the energy sweep loop (parallelized with OpenMP), and calls other modules for specific tasks.
constants_module.f90	Defines global constants and parameters for the simulation, such as physical constants, grid dimensions, time stepping details, potential parameters, and output directory names. This is the primary file to modify for changing simulation parameters.
grid_module.f90	Responsible for setting up the spatial (real-space) and momentum (k-space) grids. It defines a grid_vars_type to hold all grid-related information.
potentials_module.f90	Contains functions to calculate the values of the two diabatic potentials (E1, E2) and the real and imaginary parts of their coupling (V_real, V_imag) at any given (x,y) coordinate.
hamiltonian_module.f90	Constructs the 2x2 diabatic Hamiltonian matrix at each grid point using the functions from potentials_module.f90.
cap_module.f90	Implements the Complex Absorbing Potential (CAP) applied at the grid boundaries to prevent unphysical reflections of the wavepacket.
berry_curvature_module.f90	Calculates the Berry curvature of the lower adiabatic potential energy surface. It involves diagonalizing the diabatic Hamiltonian at multiple points (using LAPACK's ZHEEV routine) and then computing the curvature via finite differences of the eigenvectors.
fft_custom_module.f90	Provides a custom implementation of 2D Fast Fourier Transform (FFT) and its inverse (IFFT), along with a function to generate k-space frequency arrays (fftfreq equivalent).
wavepacket_init_module.f90	Initializes the starting wavepacket (Gaussian form) on the grid for a given collision energy.
propagator_module.f90	Contains the subroutines for the time propagation of the wavepacket. It implements the kinetic energy propagator (using FFTs) and the potential energy propagator (using the Crank-Nicolson scheme for the 2x2 system).
observables_module.f90	Calculates various physical observables from the wavepacket at a given time, such as total norm, expectation values of x and y, populations on diabatic states, and population-based transmission into left/right product channels and spin selectivity.
flux_module.f90	Calculates probability current densities and analyzes the flux of the wavepacket through defined planes to determine flux-based transmission probabilities and spin selectivity.
linalg_utils_module.f90	Contains small linear algebra utility functions, currently an explicit 2x2 complex matrix inversion routine used in the potential propagator.
io_module.f90	Handles input/output operations, primarily writing simulation data (adiabatic surfaces, Berry curvature, time-dependent observables, energy sweep summaries) to CSV files and managing the output directory structure.

Prerequisites and Dependencies

To compile and run this simulation, you will need:

• Fortran Compiler: The code is written in modern Fortran (F90/F95 and later features). gfortran is recommended and the provided Makefile is configured for it.
• LAPACK and BLAS Libraries: The Berry curvature calculation relies on ZHEEV from LAPACK for matrix diagonalization.
  • On Debian/Ubuntu: sudo apt-get install liblapack-dev libblas-dev
  • On Fedora/CentOS: sudo dnf install lapack-devel blas-devel
• OpenMP Support: The energy sweep loop in main_crank_prop.f90 is parallelized using OpenMP. Your Fortran compiler must support OpenMP (most modern versions of gfortran do).
• Make Utility: The make utility is used to build the project using the provided Makefile.

Compilation

A Makefile is provided to simplify the compilation process:

1. Navigate to the Source Directory: Open a terminal and change to the directory containing all the .f90 source files and the Makefile.

2. Compile: Run the make command:

   make

   This will compile all Fortran modules and the main program, then link them with the necessary libraries (LAPACK, BLAS, OpenMP) to create an executable file named crank_prop_fortran.

3. Clean Compiled Files: To remove all compiled object files (.o), module files (.mod), the executable, and the main output directory, run:

   make clean

4. Debug Build (Optional): For a build with debugging flags (e.g., for use with gdb), run:

   make debug_build

If you encounter issues with LAPACK/BLAS linking, you might need to adjust the LAPACK_LIBS variable in the Makefile to match your system's specific library paths and names.

Running the Simulation

Once the compilation is successful, you can run the simulation by executing the compiled program:

./crank_prop_fortran

The program will:

• Print status messages to the console.
• Create an output directory (default: output_fortran_crank_prop in the execution directory).
• Inside this directory, it will save:
  • CSV files for the lower and upper adiabatic potential energy surfaces.
  • A CSV file for the 2D Berry curvature map.
  • A CSV file for line cuts of the Berry curvature along x=0 and y=0.
  • For each simulated collision energy:
    • A subdirectory named energy_E.EEEEEE.
    • Inside this subdirectory, CSV files for time-dependent population-based observables (channel_data_E_...csv) and flux-based observables (flux_channel_data_E_...csv).
  • Summary CSV files in the main output directory, aggregating final transmission and selectivity values across all simulated energies.

The simulation can be computationally intensive, especially for large grids and many time steps. The OpenMP parallelization helps speed up the energy sweep.

Modifying Simulation Parameters

Most key simulation parameters are defined in constants_module.f90. To change aspects of the simulation, edit this file before recompiling. Key parameters include:

• System Parameters: A_param, omega_param, M_param, epsilon1_param, epsilon2_param, r0_param, mu_param, lambda_param (parameters for the diabatic potentials and coupling).
• Grid & Domain: nx_param, ny_param (grid points), xmin_param, xmax_param, ymin_param, ymax_param (spatial domain limits).
• Time Stepping: dt_param (time step size), num_steps_param (total number of time steps), save_every_param (frequency of saving time-dependent data).
• CAP Parameters: cap_strength_param, cap_width_x_param, cap_width_y_param.
• Output Directory: output_dir_base_name.
• Energy Sweep: num_energy_values_param, start_energy_factor, step_energy_factor (to define the range start_energy_factor*omega_param to (start_energy_factor + (num_energy_values_param-1)*step_energy_factor)*omega_param).

After modifying constants_module.f90 (or any other source file), you must recompile the project using make.

Important Note on FFT Implementation and Grid Dimensions

The current custom FFT implementation in fft_custom_module.f90 is a basic radix-2 Cooley-Tukey algorithm. This algorithm requires that the grid dimensions nx_param and ny_param (defined in constants_module.f90) be powers of 2 (e.g., 128, 256, 512, 1024, 2048, 4096, 8192).

If nx_param or ny_param are not powers of 2, the fft_1d_iterative subroutine will call STOP with an error message.

For production runs or if exact non-power-of-2 dimensions are critical:

• Consider replacing the custom FFT with a more robust library like FFTW (Fastest Fourier Transform in the West). This would involve modifying fft_custom_module.f90 to interface with FFTW and updating the Makefile to link against the FFTW library.

Output Data

The simulation generates several CSV (Comma Separated Values) files that can be analyzed and visualized using external tools like Python with Matplotlib, R, or dedicated CSV viewers.

The main outputs include:

• Adiabatic potential energy surfaces
• Berry curvature maps
• Time-dependent observables for each simulated energy
• Summary data across all energies