multidip.f90

Go to the documentation of this file.

! Copyright 2020
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-out (UKRmol+ suite).
!
!     UKRmol-out is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-out is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-out (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
!
!> \mainpage
!>
!> \brief   Calculate multi-photon dipole elements
!> \author  J Benda
!> \date    2020
!>
!> \section info Description
!>
!> Computes the multi-photon transition elements of a given order and the corresponding generalized
!> cross section. Proceeeds in the following way (see also the call graph in \ref multidip):
!>   1. The required data are read from specified K-matrix files (see \ref multidip_io::read_kmatrices)
!>      and the RMT molecular_data file (see \ref multidip_io::read_molecular_data).
!>      The K-matrix files should ideally contain the full K-matrices, i.e. including the closed
!>      channel contributions at the R-matrix radius. This can be achieved by setting `IKTYPE = 1`
!>      in R-solve. When only the (default) open-open K-matrices are used, there is a possibility
!>      of below-threshold spurious spikes in the cross sections for multi-channel calculations.
!>   2. For all intermediate orders, the intermediate states are calculated at all scattering energies
!>      (see \ref multidip_routines::solve_intermediate_state).
!>   3. The final stationary photoionization wave function is calculated at all scattering energies
!>      (see \ref multidip_routines::extract_dipole_elements).
!>   4. Transition matrix elements and photoionisation cross section are obtained
!>      (see \ref multidip_routines::calculate_pw_transition_elements).
!>
!> While the implemented method works for an arbitrary perturbation order, there are some caveats that restrict
!> the applicability of this program:
!>   1. Asymptotic series is used in place of Coulomb-Hankel functions, which does not work
!>      well close to thresholds.
!>   2. For (N + M)-photon ionisation, i.e. when M photons are absorbed in the continuum, the
!>      asymptotic complexity is proportional to O(M!). This practically restricts the tractable
!>      orders to at most 4 absorptions in the continuum.
!>
!> Multidip can either reuse existing photoionization wave function coefficient files written by RSOLVE
!> (see \ref multidip_io::read_wfncoeffs) or calculate the wave function coefficients on the fly
!> (see \ref multidip_routines::apply_Ak_coefficients_multidip).
!>
!>
!> \section cmd Command line options
!>
!> The program understands the following command line options:
!> \verbatim
!>    --help          display brief help and exit
!>    --input FILE    use FILE as the input file (instead of standard input)
!>    --test          run the internal sanity checks and exit
!> \endverbatim
!>
!> \section params Input parameters
!>
!> When the command line option `--input` is not used, MULTIDIP will read the input parameters from the
!> standard input. The input consists of the input namelist `&mdip`, which is expected to have the following form:
!>
!> \verbatim
!>   &mdip
!>     order = 2,                       ! two-photon transitions
!>
!>     rmt_data = 'molecular_data',     ! name of the molecular RMT data file (default is 'molecular_data')
!>
!>     lusct = 800,801,802,803,         ! photoionization Ak coefficient files to read (optional)
!>     lukmt = 190,191,192,194,         ! K-matrix files to read (ideally, should be produced in RSOLVE with IKTYPE = 1)
!>     nkset = 1,1,1,1,                 ! which set from the K-matrix file to use (default is to use the first set)
!>     lubnd = 500,                     ! unit with inner region expansion coefficients produced by BOUND (optional)
!>
!>     polar(:,1) = (0,0),(0,0),(1,0),  ! polarisation of the first absorbed photon (default is 0, i.e. average over polarisations)
!>     polar(:,2) = (0,0),(0,0),(1,0),  ! polarisation of the second absorbed photon (default is 0, i.e. average over polarisations)
!>
!>     omega(1) = 0,                    ! fixed energy in a.u. of the first absorbed photon or 0 for auto (default is 0)
!>     omega(2) = 0,                    ! fixed energy in a.u. of the second absorbed photon or 0 for auto (default is 0)
!>
!>     first_IP = 0,                    ! override ground state energy to yield given IP in a.u. (default is 0 = not used)
!>     rasym = 100,                     ! asymptotic radius for analytic integration of Coulomb integrals (default 100)
!>     verbose = .false.,               ! print detailed information about the execution (default is .false.)
!>     mpiio = .false.,                 ! run in HPC mode (requires MPI and ScaLAPACK, default is .false.)
!>     gpu = .false.,                   ! calculate R-matrices on GPU (requires OpenCL and CLBlast, default is .false.)
!>
!>     lab_polar = 1,1                  !lab-frame polarisations of the involved photons (default is linear, i.e. 0,0)
!>     lu_pw_dipoles = 410              !if order == 1 then the partial wave dipoles will be written to files 
!>                                      ! fort.41X, where X = MGVN, in RSOLVE format.
!>   /
!> \endverbatim
!>
!> The complete list of possible parameters and their default values are listed in documentation of the module \ref multidip_params.
!>
!> The code will use the irreducible representations specified by the K-matrix files. These need to be present also in the
!> file *molecular_data*. The user is responsible for providing all relevant K-matrix files (irreducible representations)
!> for the specified photon polarisations. Otherwise, the results may be wrong. The order of the Ak-files and K-matrices
!> has to be the same, but does not follow any specific scheme. However, the first K-matrix file given in `lukmt` defines
!> the symmetry of the initial state.
!>
!> The inner region expansion of the initial state may be provided as the output file of BOUND by means of the keyword `lubnd`.
!> This expansion is then automatically continued to the outer region by MULTIDIP and taken in account in further calculations.
!> This increases the computational complexity in the same way as raising the amount of absorbed photons by one. If no bound
!> state coefficients are provided, the program will simply use the lower-energy inner region eigenstate in the symmetry specified
!> by the first K-matrix file.
!>
!> The polarisation of \f$ k \f$-th photon is given as three complex-valued components of a unit Cartesian vector:
!> \f$ (\mathrm{Re}\,\epsilon_x^{(k)}, \mathrm{Im}\,\epsilon_x^{(k)}),
!>     (\mathrm{Re}\,\epsilon_y^{(k)}, \mathrm{Im}\,\epsilon_y^{(k)}),
!>     (\mathrm{Re}\,\epsilon_z^{(k)}, \mathrm{Im}\,\epsilon_z^{(k)}) \f$.
!> When a polarisation vector has only real components (as in the above example), some compilers (like ifort, but not gfortran)
!> will allow omission of the parentheses and specifying just the real part. When polarisation of some photon is specified as zero,
!> no partial wave multipole moments will be written to files at the end of the calculation.
!>
!> Some fixed absorbed photon energies can be specified on input using the `omega` parameter. There must always be at least one
!> photon with variable energy (corresponds to `omega = 0`). The energies of the variable-energy photons will be then equally
!> assigned so that the total energy of the system corresponds to the total energies on input (in the K-matrix files). By default
!> `omega = 0` for all photons, so all photons have variable energy and will correspond to the same wave length.
!>
!> If non-zero, the parameter `rasym` defines the asymptotic radius. It has to be larger than the radius of the inner region.
!> Radial integrals of Coulomb functions are integrated numerically between the inner region boundary and this radius,
!> while asymptotic expansion and analytic integration is used from the asymptotic radius to infinity. The value of 100 atomic
!> units is typically a good choice. Note that numerical integration is implemented only for two-photon processes and this
!> parameter will be presently ignored for higher-order processes.
!>
!> \section outputs Output files
!>
!> When the polarisations of the photons are given, the program produces the file *gen_photo_xsec.txt* with the generalized cross
!> section for oriented molecule. The first column is the energy of the first photon in eV, the second column is the total cross
!> section in atomic units (area^N time^(N-1)), further columns are partial cross sections associated with transition to a
!> specific final ionic state. See \ref multidip_routines::calculate_pw_transition_elements for details.
!>
!> Likewise, if all photon polarisations are given, the full partial wave N-photon transition matrix elements will be written to
!> files *pwdip-M-(I,L,M).txt* and *pwdip-M-(I,L,M)+cphase.txt*, where M is the irreducible representation (0-based), I is
!> the (1-based) index of the residual ion, and L and M are the angular momentum quantum numbers of the emitted electron's partial
!> wave. The file with "+cphase" ending contains the additional partial wave phase factors
!> \f$ \mathrm{i}^\ell \exp(-\mathrm{i} \sigma_\ell) \f$. The first column
!> of these files is the real part of the matrix element, the second is the imaginary part. The total energies are not present,
!> but correspond to the same energies as given by the input K-matrix files.
!>
!> Irrespective of the polarisations in the input namelist, MULTIDIP will produce the files *gen_<n>photo_beta_<L>.txt* with the
!> parameters \f$ \beta_L \f$ of the molecular-orientatio-averaged photoelectron angular distribution,
!> \f[
!>     \frac{\mathrm{d}\sigma^{(p_1\dots p_n)}}{\mathrm{d}\Omega}
!>     = \frac{1}{4\pi} \beta_0 \left(1 + \beta_1 P_1(\cos\theta) + \dots + \beta_{2n} P_{2n}(\cos\theta) \right) ,
!> \f]
!> where \f$ \beta_0 = \sigma^{(p_1\dots p_n)} \f$ is the integrated partial cross section and \f$ \beta_{L > 0} \f$
!> is the dimensionless asymmetry parameter for angular momentum \f$ L \f$. Currently, the laboratory-frame polarisations
!> \f$ p_1,\dots,p_n \f$ are all hardcoded to zero (linear polarisation). Files that would contain zeros only are not written.
!> The layout of these files is the same as of *gen_photo_xsec.txt*.
!>
!> Finally, if some fully custom processing of the raw transition dipoles is required, it is possible to let them MULTIDIP write all
!> to disk by use of the keyword `raw`. When set to the string "xyz" or "sph" or "both", raw transition dipoles in the selected
!> angular basis will be written to files *rawdip-CHAIN-(I,L,M).txt*. Here "CHAIN" will be the absorption history of the transition,
!> for instance "xx" if that raw dipole corresponds to absorption of two x-polarized photons in the Cartesian basis, or "mp0" if it
!> corresponds to absorption of \f$ m = 0 \f$, \f$ m = +1 \f$ and \f$ m = -1 \f$ photons in the spherical basis. The first absorbed
!> photon corresponds to the right-most index, as is customary in the typical bra-ket notation. The files contain
!> real and imaginary parts of the dipoles for all scattering energies. The utility programs \ref multidip_smooth_rawdips and
!> \ref multidip_rabitt_delays work with these files.
!>
!>
!> \section opts Configuration options
!>
!> | CMake option          | default value   |
!> | ---                   | ---             |
!> | `WITH_CLBLAST`        | `OFF`           |
!> | `WITH_GSL`            | `OFF`           |
!> | `WITH_MMAP`           | `ON`            |
!> | `SCALAPACK_LIBRARIES` | (empty)         |
!>
!> This program can optionally use special functions (the complex gamma function, Coulomb functions and the confluent
!> hypergeometric functions) from the GNU Scientific Library instead of the UKRmol-out library. For that, pass the
!> option `-D WITH_GSL=ON` to CMake during configuration.
!>
!> When configured with the `-D WITH_MMAP=ON` option, the inner region dipole matrices will not be read
!> into memory, but only mapped into virtual memory. (This requires presence of such functionality
!> in GBTOlib.)
!>
!> This program benefits from optimized and/or threaded BLAS. Some sections of the code are also explicitly threaded (using
!> OpenMP). Moreover, when compiled with the MPI support (this is controlled by GBTOlib level), the loop over energies
!> (as well as some other pieces of the code) will be distributed among the parallel images in the round-robin fashion.
!>
!> When the support for parallel MPI is enabled during compilation, it is also possible to enable use of MPI-IO and ScaLAPACK.
!> MPI-IO is used to read the large inner region dipole matrices. ScaLAPACK is used to redistribute solutions and to multiply them
!> with the inner dipoles. This MPI-IO/ScaLAPACK support will be compiled in if the option `SCALAPACK_LIBRARIES` is passed to
!> CMake, and the code will be activated by setting `mpiio = .true.` in the program's input namelist.
!>
!> The computation of R-matrices (R = w^T [E - H]^(-1) w), which is dominated by the multiplication of large constant real boundary
!> amplitude matrices to produce a much smaller output matrix, can be offloaded to a GPU. This requires compiling UKRmol-out with
!> the library CLBlast by specifying `WITH_CLBLAST=ON` and also `CLBLAST_INCLUDE_DIRS` and `CLBLAS_LIBRARIES` when running CMake.
!> The actual GPU mode is then enabled by specifying `gpu = .true.` in the input namelist. In serial mode the OpenCL platform and
!> OpenCL device index to use (both default 0) can be set by the environment variables `OCL_PLATFORM` and `OCL_DEVICE`. In parallel
!> mode the latter variable is ignored and the parallel processes are associated with the available OpenCL devices within the
!> selected platform in the round-robin way.
!>
program multidip
 
    use multidip_routines, only: multidip_main
 
    implicit none
 
    call multidip_main
 
end program multidip