! This Source Code Form is subject to the terms of the Mozilla Public
! License, v. 2.0. If a copy of the MPL was not distributed with this
! file, You can obtain one at http://mozilla.org/MPL/2.0/.
module mbd
!! High-level Fortran API.

use mbd_constants
use mbd_defaults
use mbd_version
use mbd_damping, only: damping_t
use mbd_formulas, only: scale_with_ratio
use mbd_geom, only: geom_t
use mbd_gradients, only: grad_request_t, grad_t
use mbd_methods, only: get_mbd_energy, get_mbd_scs_energy
use mbd_ts, only: get_ts_energy
use mbd_utils, only: result_t, exception_t, printer_i
use mbd_vdw_param, only: ts_vdw_params, tssurf_vdw_params, species_index

implicit none

private

public :: MBD_VERSION_MAJOR, MBD_VERSION_MINOR, MBD_VERSION_PATCH, MBD_VERSION_SUFFIX
public :: MBD_EXC_NEG_EIGVALS, MBD_EXC_NEG_POL, MBD_EXC_LINALG, MBD_EXC_UNIMPL, &
    MBD_EXC_DAMPING, MBD_EXC_INPUT

type, public :: mbd_input_t
    !! Contains user input to an MBD calculation.
    character(len=30) :: method = 'mbd-rsscs'
        !! VdW method to use to calculate energy and gradients.
        !!
        !! - `mbd-rsscs`: The MBD@rsSCS method.
        !! - `mbd-nl`: The MBD-NL method.
        !! - `ts`: The TS method.
        !! - `mbd`: Generic MBD method (without any screening).
    integer :: comm = -1
        !! MPI communicator, as a plain integer handle.
        !!
        !! Only used when compiled with MPI. Leave as is to use the
        !! MPI_COMM_WORLD communicator. The integer handle is accepted
        !! regardless of whether the library was built against the
        !! `mpi_f08` module; callers using `mpi_f08` should pass the
        !! integer handle of their communicator (e.g. `mycomm%mpi_val`).
    integer :: max_atoms_per_block = MAX_ATOMS_PER_BLOCK
        !! Number of atoms per block in a BLACS grid.
    integer :: log_level = MBD_LOG_LVL_INFO
        !! Level of printing
    procedure(printer_i), nopass, pointer :: printer => null()
        !! If assigned, will be used for logging
    logical :: calculate_forces = .true.
        !! Whether to calculate forces.
    logical :: calculate_vdw_params_gradients = .false.
        !! Whether to calculate gradients of energy w.r.t. vdW parameters
    logical :: calculate_spectrum = .false.
        !! Whether to keep MBD eigenvalues.
    logical :: do_rpa = .false.
        !! Whether to evalulate the MBD energy as an RPA integral over frequency.
    logical :: rpa_orders = .false.
        !! Whether to calculate individual RPA orders
    logical :: rpa_rescale_eigs = .false.
        !! Whether to rescale RPA eigenvalues as in 10.1021/acs.jctc.6b00925.
    integer :: n_omega_grid = N_FREQUENCY_GRID
        !! Number of imaginary frequency grid points.
    real(dp) :: k_grid_shift = K_GRID_SHIFT
        !! Off-\(\Gamma\) shift of the \(k\)-point grid in units of
        !! inter-\(k\)-point distance.
    logical :: zero_negative_eigvals = .false.
        !! Whether to zero out negative eigenvalues.
    character(len=20) :: xc = ''
        !! XC functional for automatic setting of damping parameters.
    real(dp) :: ts_d = TS_DAMPING_D
        !! TS damping parameter \(d\).
    real(dp) :: ts_sr = -1
        !! Custom TS damping parameter \(s_R\).
        !!
        !! Leave as is to use a value based on the XC functional.
    real(dp) :: mbd_a = MBD_DAMPING_A
        !! MBD damping parameter \(a\).
    real(dp) :: mbd_beta = -1
        !! Custom MBD damping parameter \(\beta\).
        !!
        !! Leave as is to use a value based on the XC functional.
    character(len=10) :: vdw_params_kind = 'ts'
        !! Which free-atom reference vdW parameters to use for scaling.
        !!
        !! - `ts`: Values from original TS method.
        !! - `tssurf`: Values from the TS\(^\text{surf}\) approach.
    character(len=3), allocatable :: atom_types(:)
        !! (\(N\)) Atom types used for picking free-atom reference values.
    real(dp), allocatable :: free_values(:, :)
        !! (\(N\times3\), a.u.) Custom free-atom vdW paramters to use for
        !! scaling.
        !!
        !! Columns contain static polarizabilities, C6 coefficients, and vdW
        !! radii.
    real(dp), allocatable :: coords(:, :)
        !! (\(3\times N\), a.u.) Atomic coordinates.
    real(dp), allocatable :: lattice_vectors(:, :)
        !! (\(3\times 3\), a.u.) Lattice vectors in columns, unallocated if not
        !! periodic.
    integer :: k_grid(3) = [-1, -1, -1]
        !! Number of \(k\)-points along reciprocal axes.
    character(len=10) :: parallel_mode = 'auto'
        !! Parallelization scheme.
        !!
        !! - `auto`: Pick based on system system size and number of \(k\)-points.
        !! - `kpoints`: Parallelize over \(k\)-points.
        !! - `atoms`: Parallelize over atom pairs.
end type

type, public :: mbd_calc_t
    !! Represents an MBD calculation.
    private
    type(geom_t) :: geom
    type(damping_t) :: damp
    real(dp), allocatable :: alpha_0(:)
    real(dp), allocatable :: C6(:)
    character(len=30) :: method
    type(result_t) :: results
    type(grad_t) :: dalpha_0, dC6, dr_vdw
    logical :: calculate_gradients
    logical :: calculate_vdw_params_gradients
    real(dp), allocatable :: free_values(:, :)
    character(len=30) :: vdw_params_update
contains
    procedure :: init => mbd_calc_init
    procedure :: destroy => mbd_calc_destroy
    procedure :: switch_forces => mbd_calc_switch_forces
    procedure :: update_coords => mbd_calc_update_coords
    procedure :: update_lattice_vectors => mbd_calc_update_lattice_vectors
    procedure :: update_vdw_params_custom => mbd_calc_update_vdw_params_custom
    procedure :: update_vdw_params_from_ratios => mbd_calc_update_vdw_params_from_ratios
    procedure :: update_vdw_params_nl => mbd_calc_update_vdw_params_nl
    procedure :: evaluate_vdw_method => mbd_calc_evaluate_vdw_method
    procedure :: get_gradients => mbd_calc_get_gradients
    procedure :: get_vdw_params_ratios_gradients => mbd_calc_get_vdw_params_ratios_gradients
    procedure :: get_lattice_derivs => mbd_calc_get_lattice_derivs
    procedure :: get_lattice_stress => mbd_calc_get_lattice_stress
    procedure :: get_spectrum_modes => mbd_calc_get_spectrum_modes
    procedure :: get_rpa_orders => mbd_calc_get_rpa_orders
    procedure :: get_exception => mbd_calc_get_exception
end type

contains

subroutine mbd_calc_init(this, input)
    !! Initialize an MBD calculation from an MBD input.
    class(mbd_calc_t), target, intent(inout) :: this
    type(mbd_input_t), intent(in) :: input
        !! MBD input.

#ifdef WITH_MPI
    if (input%comm /= -1) then
#   ifdef WITH_MPIF08
        this%geom%mpi_comm%mpi_val = input%comm
#   else
        this%geom%mpi_comm = input%comm
#   endif
    end if
#endif
#ifdef WITH_SCALAPACK
    this%geom%max_atoms_per_block = input%max_atoms_per_block
#endif
    this%method = input%method
    this%calculate_gradients = input%calculate_forces
    this%calculate_vdw_params_gradients = input%calculate_vdw_params_gradients
    this%geom%get_eigs = input%calculate_spectrum
    this%geom%get_modes = input%calculate_spectrum
    this%geom%do_rpa = input%do_rpa
    this%geom%get_rpa_orders = input%rpa_orders
    this%geom%param%rpa_rescale_eigs = input%rpa_rescale_eigs
    this%geom%param%n_freq = input%n_omega_grid
    this%geom%param%k_grid_shift = input%k_grid_shift
    this%geom%param%zero_negative_eigvals = input%zero_negative_eigvals
    if (.not. all(input%k_grid == -1)) this%geom%k_grid = input%k_grid
    this%geom%coords = input%coords
    if (allocated(input%lattice_vectors)) then
        if (input%method /= 'ts' .and. .not. allocated(this%geom%k_grid)) then
            this%geom%exc = exception_t( &
                MBD_EXC_INPUT, &
                'calc%init()', &
                'Lattice vectors present but no k-grid specified' &
            )
            return
        end if
        this%geom%lattice = input%lattice_vectors
    end if
    this%geom%parallel_mode = input%parallel_mode
    if (associated(input%printer)) this%geom%log%printer => input%printer
    this%geom%log%level = input%log_level
    call this%geom%init()
    if (allocated(input%free_values)) then
        this%free_values = input%free_values
    else
        select case (input%vdw_params_kind)
        case ('ts')
            this%free_values = ts_vdw_params(:, species_index(input%atom_types))
        case ('tssurf')
            this%free_values = tssurf_vdw_params(:, species_index(input%atom_types))
        end select
    end if
    if (input%xc == '') then
        this%damp%beta = input%mbd_beta
        this%damp%a = input%mbd_a
        this%damp%ts_d = input%ts_d
        this%damp%ts_sr = input%ts_sr
        select case (input%method)
        case ('ts')
            if (input%ts_sr < 0) then
                this%geom%exc%code = MBD_EXC_DAMPING
                this%geom%exc%msg = 'Damping parameter S_r for TS not specified'
            end if
        case default
            if (input%mbd_beta < 0) then
                this%geom%exc%code = MBD_EXC_DAMPING
                this%geom%exc%msg = 'Damping parameter beta for MBD not specified'
            end if
        end select
    else
        this%geom%exc = this%damp%set_params_from_xc(input%xc, input%method)
    end if
    if (this%geom%has_exc()) return
end subroutine

subroutine mbd_calc_destroy(this)
    !! Finalize an MBD calculation.
    class(mbd_calc_t), target, intent(inout) :: this

    call this%geom%destroy()
end subroutine

subroutine mbd_calc_switch_forces(this, forces)
    !! Update whether to calculate forces.
    class(mbd_calc_t), intent(inout) :: this
    logical, intent(in) :: forces
        !! Whether to calcualte forces.

    this%calculate_gradients = forces
end subroutine

subroutine mbd_calc_update_coords(this, coords)
    !! Update atomic coordinates.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(in) :: coords(:, :)
        !! (\(3\times N\), a.u.) New atomic coordinates.

    this%geom%coords = coords
end subroutine

subroutine mbd_calc_update_lattice_vectors(this, latt_vecs)
    !! Update unit-cell lattice vectors.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(in) :: latt_vecs(:, :)
        !! (\(3\times 3\), a.u.) New lattice vectors in columns.

    this%geom%lattice = latt_vecs
end subroutine

subroutine mbd_calc_update_vdw_params_custom(this, alpha_0, C6, r_vdw)
    !! Update vdW parameters in a custom way.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(in) :: alpha_0(:)
        !! (a.u.) New atomic static polarizabilities.
    real(dp), intent(in) :: C6(:)
        !! (a.u.) New atomic \(C_6\) coefficients.
    real(dp), intent(in) :: r_vdw(:)
        !! (a.u.) New atomic vdW radii.

    this%alpha_0 = alpha_0
    this%C6 = C6
    this%damp%r_vdw = r_vdw
    this%vdw_params_update = 'custom'
end subroutine

subroutine mbd_calc_update_vdw_params_from_ratios(this, ratios)
    !! Update vdW parameters based on scaling of free-atom values.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(in) :: ratios(:)
        !! Ratios of atomic volumes in the system and in vacuum.

    real(dp), allocatable :: ones(:)
    type(grad_request_t) :: grad

    allocate (ones(size(ratios)), source=1d0)
    grad%dV = this%calculate_vdw_params_gradients
    ! Alias free_values before slicing it into the scale_with_ratio calls.
    ! Passing an array section of an allocatable rank-2 component of the
    ! polymorphic `this` straight into an array-valued function makes gfortran 16
    ! emit an -fcheck=bounds check that reads the component descriptor through a
    ! null pointer -- a bogus bounds error at -O0, a segfault at -O1/-Og
    ! (libmbd/libmbd#127). Associating the component evaluates its descriptor
    ! once and sidesteps the bad check; it is correct and idiomatic on every
    ! compiler, so no version guard is needed. Do not fold the alias back into
    ! `this%free_values(i, :)` while gfortran 16 is supported.
    !
    ! Scope of the bug: reproduced on the gfortran 16 series only -- 16.0 trunk
    ! (r16-8100) and the 16.1.0 release (both x86-64 and aarch64). gfortran 13,
    ! 14 and 15 are unaffected, and it is already fixed on the 17.0 trunk, so the
    ! alias can be removed once gfortran 16 is no longer supported. No upstream
    ! GCC bug was filed (it is fixed on trunk); a minimal standalone reproducer
    ! is kept in the PR discussion for libmbd/libmbd#127.
    associate (free_values => this%free_values)
    this%alpha_0 = scale_with_ratio( &
        free_values(1, :), ratios, ones, 1d0, this%dalpha_0, grad &
    )
    this%C6 = scale_with_ratio( &
        free_values(2, :), ratios, ones, 2d0, this%dC6, grad &
    )
    this%damp%r_vdw = scale_with_ratio( &
        free_values(3, :), ratios, ones, 1d0 / 3, this%dr_vdw, grad &
    )
    end associate
    this%vdw_params_update = 'ratios'
end subroutine

subroutine mbd_calc_get_vdw_params_ratios_gradients(this, dE_dratios)
    !! Get gradients of the energy w.r.t. Hirshfeld ratios if they were
    !! requested in the MBD input.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(out) :: dE_dratios(:)
        !! Gradients of the energy w.r.t. Hirshfeld ratios.

    if (this%vdw_params_update /= 'ratios') return
    dE_dratios = ( &
        this%results%dE%dalpha * this%dalpha_0%dV &
        + this%results%dE%dC6 * this%dC6%dV &
        + this%results%dE%dr_vdw * this%dr_vdw%dV &
     )
end subroutine

subroutine mbd_calc_update_vdw_params_nl(this, alpha_0_ratios, C6_ratios)
    !! Update vdW parameters for the MBD-NL method.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(in) :: alpha_0_ratios(:)
        !! Ratios of free-atom exact static polarizabilities and those from the
        !! VV functional.
    real(dp), intent(in) :: C6_ratios(:)
        !! Ratios of free-atom exact \(C_6\) coefficients and those from the VV
        !! functional.

    this%alpha_0 = this%free_values(1, :) * alpha_0_ratios
    this%C6 = this%free_values(2, :) * C6_ratios
    this%damp%r_vdw = 2.5d0 * this%free_values(1, :)**(1d0 / 7) * alpha_0_ratios**(1d0 / 3)
    this%vdw_params_update = 'nl'
end subroutine

subroutine mbd_calc_evaluate_vdw_method(this, energy)
    !! Evaluate a given vdW method for a given system and vdW parameters,
    !! retrieve energy.
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(out) :: energy
        !! (a.u.) VdW energy.

    type(grad_request_t) :: grad

    if (this%calculate_gradients) then
        grad%dcoords = .true.
        if (allocated(this%geom%lattice)) grad%dlattice = .true.
    end if
    if (this%calculate_vdw_params_gradients) then
        grad%dalpha = .true.
        grad%dC6 = .true.
        grad%dr_vdw = .true.
    end if
    select case (this%method)
    case ('mbd', 'mbd-nl')
        this%damp%version = 'fermi,dip'
        this%results = get_mbd_energy( &
            this%geom, this%alpha_0, this%C6, this%damp, grad &
        )
        energy = this%results%energy
    case ('mbd-rsscs')
        this%results = get_mbd_scs_energy( &
            this%geom, 'rsscs', this%alpha_0, this%C6, this%damp, grad &
        )
        energy = this%results%energy
    case ('ts')
        this%damp%version = 'fermi'
        this%results = get_ts_energy( &
            this%geom, this%alpha_0, this%C6, this%damp, grad &
        )
        energy = this%results%energy
    end select
    if (this%geom%log%level <= MBD_LOG_LVL_DEBUG) call this%geom%timer%print()
end subroutine

subroutine mbd_calc_get_gradients(this, gradients)  ! 3 by N  dE/dR
    !! Retrieve nuclear energy gradients if they were requested in the MBD
    !! input.
    !!
    !! The gradients are calculated together with the energy, so a call to this
    !! method must be preceeded by a call to
    !! [[mbd_calc_t:evaluate_vdw_method]].  For the same reason, the
    !! gradients must be requested prior to this called via
    !! [[mbd_input_t:calculate_forces]].
    class(mbd_calc_t), intent(in) :: this
    real(dp), intent(out) :: gradients(:, :)
        !! (\(3\times N\), a.u.) Energy gradients, \(\mathrm dE/\mathrm d\mathbf
        !! R_i\), index \(i\) runs over columns.

    gradients = transpose(this%results%dE%dcoords)
end subroutine

subroutine mbd_calc_get_lattice_derivs(this, latt_derivs)
    !! Provide lattice-vector energy gradients if they were requested in the MBD
    !! input.
    !!
    !! The gradients are actually calculated together with the energy, so a call
    !! to this method must be preceeded by a call to
    !! [[mbd_calc_t:evaluate_vdw_method]].  For the same reason, the
    !! gradients must be requested prior to this called via
    !! [[mbd_input_t:calculate_forces]].
    class(mbd_calc_t), intent(in) :: this
    real(dp), intent(out) :: latt_derivs(:, :)
        !! (\(3\times 3\), a.u.) Energy gradients, \(\mathrm dE/\mathrm d\mathbf
        !! a_i\), index \(i\) runs over columns.

    latt_derivs = transpose(this%results%dE%dlattice)
end subroutine

subroutine mbd_calc_get_lattice_stress(this, stress)
    !! Provide stress tensor of the lattice.
    !!
    !! This is a utility function wrapping [[mbd_calc_t:get_lattice_derivs]].
    !! The lattice vector gradients are coverted to the stress tensor.
    class(mbd_calc_t), intent(in) :: this
    real(dp), intent(out) :: stress(:, :)
        !! (\(3\times 3\), a.u.) Stress tensor.

    stress = ( &
        matmul(this%geom%lattice, this%results%dE%dlattice) &
        + matmul(this%geom%coords, this%results%dE%dcoords) &
    )
end subroutine

subroutine mbd_calc_get_spectrum_modes(this, spectrum, modes)
    !! Provide MBD spectrum if it was requested in the MBD input.
    !!
    !! The spectrum is actually calculated together with the energy, so a call
    !! to this method must be preceeded by a call to
    !! [[mbd_calc_t:evaluate_vdw_method]].  For the same reason, the
    !! spectrum must be requested prior to this call via
    !! [[mbd_input_t:calculate_spectrum]].
    class(mbd_calc_t), intent(inout) :: this
    real(dp), intent(out) :: spectrum(:)
        !! (\(3N\), a.u.) Energies (frequencies) of coupled MBD modues,
        !! \(\omega_i\).
    real(dp), intent(out), allocatable, optional :: modes(:, :)
        !! (\(3N\times 3N\)) Coupled-mode wave functions (MBD eigenstates),
        !! \(\psi_j\), in the basis of uncoupled states,
        !! \(C_{ij}=\langle\phi_i|\psi_j\rangle\), index \(j\) runs over
        !! columns.
        !!
        !! To save memory, the argument must be allocatable, and the method
        !! transfers allocation from the internal state to the argument. For
        !! this reason, the method can be called only once wih this optional
        !! argument per calculation.

    spectrum = this%results%mode_eigs
    if (present(modes)) call move_alloc(this%results%modes, modes)
end subroutine

subroutine mbd_calc_get_rpa_orders(this, rpa_orders)
    !! Provide RPA orders if they were requested in the MBD input.
    !!
    !! The orders are actually calculated together with the energy, so a call
    !! to this method must be preceeded by a call to
    !! [[mbd_calc_t:evaluate_vdw_method]]. For the same reason, the
    !! spectrum must be requested prior to this call via
    !! [[mbd_input_t:do_rpa]] and [[mbd_input_t:rpa_orders]].
    class(mbd_calc_t), intent(inout) :: this
    real(dp), allocatable, intent(out) :: rpa_orders(:)
        !! (a.u.) MBD energy decomposed to RPA orders.

    rpa_orders = this%results%rpa_orders
end subroutine

subroutine mbd_calc_get_exception(this, code, origin, msg)
    !! Retrieve an exception in the MBD calculation if it occured.
    class(mbd_calc_t), intent(inout) :: this
    integer, intent(out) :: code
        !! Exception code, values defined in [[mbd_constants]].
    character(*), intent(out) :: origin
        !! Exception origin.
    character(*), intent(out) :: msg
        !! Exception message.

    code = this%geom%exc%code
    if (code == 0) return
    origin = this%geom%exc%origin
    msg = this%geom%exc%msg
    this%geom%exc%code = 0
    this%geom%exc%origin = ''
    this%geom%exc%msg = ''
end subroutine

end module
