mod_solvers Module

Module that contains the routines used to solve the polarization linear system . Currently three methods are implemented:
1. matrix inversion;
2. (preconditioned) conjugate gradients - since polarization equations are symmetric and positive definite, this is the optimal choice;
3. jacobi iterations accelerated with Pulay's direct inversion in the iterative subspace (DIIS): this is a pretty robust solver that can be use for general systems and that is less sensitive to small errors in the symmetry of the matrix.

Iterative solvers need two additional routines to be passed as arguments, namely matvec that computes a generic product and precond that computes , where is a precontioner


Uses

  • module~~mod_solvers~~UsesGraph module~mod_solvers mod_solvers module~mod_io mod_io module~mod_solvers->module~mod_io module~mod_memory mod_memory module~mod_solvers->module~mod_memory module~mod_constants mod_constants module~mod_solvers->module~mod_constants module~mod_electrostatics mod_electrostatics module~mod_solvers->module~mod_electrostatics module~mod_io->module~mod_constants module~mod_memory->module~mod_io module~mod_memory->module~mod_constants iso_c_binding iso_c_binding module~mod_memory->iso_c_binding module~mod_constants->iso_c_binding module~mod_electrostatics->module~mod_io module~mod_electrostatics->module~mod_memory module~mod_electrostatics->module~mod_constants module~mod_profiling mod_profiling module~mod_electrostatics->module~mod_profiling module~mod_adjacency_mat mod_adjacency_mat module~mod_electrostatics->module~mod_adjacency_mat module~mod_topology mod_topology module~mod_electrostatics->module~mod_topology module~mod_profiling->module~mod_io module~mod_profiling->module~mod_memory module~mod_profiling->module~mod_constants module~mod_adjacency_mat->module~mod_memory module~mod_topology->module~mod_memory module~mod_topology->module~mod_adjacency_mat

Used by

  • module~~mod_solvers~~UsedByGraph module~mod_solvers mod_solvers proc~polarization polarization proc~polarization->module~mod_solvers

Contents

Variables

Subroutines


Variables

Type Visibility Attributes Name Initial
real(kind=rp), private, parameter :: OMMP_DEFAULT_SOLVER_TOL = 1e-8_rp

Default tolerance for iterative solvers
integer(kind=ip), private, parameter :: OMMP_DEFAULT_SOLVER_ITER = 200

Default maximum number of iteration for iterative solvers
integer(kind=ip), private, parameter :: OMMP_DEFAULT_DIIS_MAX_POINTS = 20

Default maximum number of points in DIIS extrapolation

Subroutines

public subroutine inversion_solver(n, rhs, x, tmat)

Solve the linear system directly inverting the matrix: This is highly unefficient and should only be used for testing other methods of solution.

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n

Size of the matrix
real(kind=rp), intent(in), dimension(n) :: rhs

Right hand side of the linear system
real(kind=rp), intent(out), dimension(n) :: x

In output the solution of the linear system
real(kind=rp), intent(in), dimension(n, n) :: tmat

Polarization matrix TODO

public subroutine conjugate_gradient_solver(n, rhs, x, eel, matvec, precnd, arg_tol, arg_n_iter)

Conjugate gradient solver (TODO)

Read more…

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n

Size of the matrix
real(kind=rp), intent(in), dimension(n) :: rhs

Right hand side of the linear system
real(kind=rp), intent(inout), dimension(n) :: x

In input, initial guess for the solver, in output the solution
type(ommp_electrostatics_type), intent(in) :: eel

Electrostatics data structure
integer :: matvec
real :: precnd
real(kind=rp), intent(in), optional :: arg_tol

Optional convergence criterion in input, if not present OMMP_DEFAULT_SOLVER_TOL is used.
integer(kind=ip), intent(in), optional :: arg_n_iter

Optional maximum number of iterations for the solver, if not present OMMP_DEFAULT_SOLVER_ITER is used.

public subroutine jacobi_diis_solver(n, rhs, x, eel, matvec, inv_diag, arg_tol, arg_n_iter, arg_diis_max)

Routine to perform matrix-vector product

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n

Size of the matrix
real(kind=rp), intent(in), dimension(n) :: rhs

Right hand side of the linear system
real(kind=rp), intent(inout), dimension(n) :: x

In input, initial guess for the solver, in output the solution
type(ommp_electrostatics_type), intent(in) :: eel

Electrostatics data structure
integer :: matvec
real(kind=rp), intent(in), dimension(n) :: inv_diag

Element-wise inverse of diagonal of LHS matrix
real(kind=rp), intent(in), optional :: arg_tol

Optional convergence criterion in input, if not present OMMP_DEFAULT_SOLVER_TOL is used.
integer(kind=ip), intent(in), optional :: arg_n_iter

Optional maximum number of iterations for the solver, if not present OMMP_DEFAULT_SOLVER_ITER is used.
integer(kind=ip), intent(in), optional :: arg_diis_max

Optional maximum number of points for diis extrapolation, if not present OMMP_DEFAULT_DIIS_MAX_POINTS is used.

private subroutine diis(n, nmat, ndiis, x, e, b, xnew)

perform Pulay's direct inversion in the iterative subspace extrapolation:

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n
integer(kind=ip), intent(inout) :: nmat
integer(kind=ip), intent(in) :: ndiis
real(kind=rp), intent(inout), dimension(n, ndiis) :: x
real(kind=rp), intent(inout), dimension(n, ndiis) :: e
real(kind=rp), intent(inout), dimension(ndiis+1, ndiis+1) :: b
real(kind=rp), intent(inout), dimension(n) :: xnew

private subroutine makeb(n, nmat, ndiis, e, b)

assemble the DIIS B matrix:

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n
integer(kind=ip), intent(in) :: nmat
integer(kind=ip), intent(in) :: ndiis
real(kind=rp), intent(in), dimension(n, ndiis) :: e
real(kind=rp), intent(inout), dimension(ndiis+1, ndiis+1) :: b

private subroutine rmsvec(n, v, vrms, vmax)

compute root-mean-square and max norms of a vector.

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n
real(kind=rp), intent(in), dimension(n) :: v
real(kind=rp), intent(inout) :: vrms
real(kind=rp), intent(inout) :: vmax