Main driver for the calculation of induced dipoles. Takes electric field at induced dipole sites as input and -- if solver converges -- provides induced dipoles as output. Since AMOEBA requires the calculations of two sets of induced dipoles generated from two different electric fields (normally called direct (D) and polarization (P)) both electric field and induced dipoles are shaped with an extra dimension and this routine calls the solver twice to solve the two linear systems in the case of AMOEBA FF. Direct electric field and induced dipoles are stored in e(:,:,1)/ipds(:,:,1) while polarization field/dipole are stored in e(:,:,2)/ipds(:,:,2).
TODO Maybe check convergence...
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_system), | intent(inout), | target | :: | sys_obj |
Fundamental data structure for OMMP system |
|
| real(kind=rp), | intent(in), | dimension(3, sys_obj%eel%pol_atoms, sys_obj%eel%n_ipd) | :: | e |
Total electric field that induces the dipoles |
|
| integer(kind=ip), | intent(in), | optional | :: | arg_solver |
Flag for the solver to be used; optional, should be one OMMP_SOLVER_ if not provided ommp_solver_default is used. |
|
| integer(kind=ip), | intent(in), | optional | :: | arg_mvmethod |
Flag for the matrix-vector method to be used; optional, should be one of OMMP_MATV_ if not provided ommp_matv_default is used. |
|
| logical, | intent(in), | optional | :: | arg_ipd_mask(sys_obj%eel%n_ipd) |
Logical mask to skip calculation of one of the two set of dipoles in AMOEBA calculations (eg. when MM part's field is not taken into account, both P and D field are just external field, so there is no reason to compute it twice). If n_ipd == 1 this is always considered true. |
subroutine polarization(sys_obj, e, &
& arg_solver, arg_mvmethod, arg_ipd_mask)
!! Main driver for the calculation of induced dipoles.
!! Takes electric field at induced dipole sites as input and -- if
!! solver converges -- provides induced dipoles as output.
!! Since AMOEBA requires the calculations of two sets of induced dipoles
!! generated from two different electric fields (normally called direct (D)
!! and polarization (P)) both electric field and induced dipoles are shaped
!! with an extra dimension and this routine calls the solver twice to
!! solve the two linear systems in the case of AMOEBA FF. Direct electric
!! field and induced dipoles are stored in e(:,:,1)/ipds(:,:,1) while
!! polarization field/dipole are stored in e(:,:,2)/ipds(:,:,2).
use mod_solvers, only: jacobi_diis_solver, conjugate_gradient_solver, &
inversion_solver
use mod_memory, only: ip, rp, mallocate, mfree
use mod_io, only: print_matrix
use mod_profiling, only: time_pull, time_push
use mod_constants, only: OMMP_MATV_DIRECT, &
OMMP_MATV_INCORE, &
OMMP_MATV_NONE, &
OMMP_SOLVER_CG, &
OMMP_SOLVER_DIIS, &
OMMP_SOLVER_INVERSION, &
OMMP_SOLVER_NONE, &
OMMP_VERBOSE_DEBUG, &
OMMP_VERBOSE_HIGH
implicit none
type(ommp_system), target, intent(inout) :: sys_obj
!! Fundamental data structure for OMMP system
real(rp), dimension(3, sys_obj%eel%pol_atoms, sys_obj%eel%n_ipd), &
& intent(in) :: e
!! Total electric field that induces the dipoles
integer(ip), intent(in), optional :: arg_solver
!! Flag for the solver to be used; optional, should be one OMMP_SOLVER_
!! if not provided [[mod_constants:OMMP_SOLVER_DEFAULT]] is used.
integer(ip), intent(in), optional :: arg_mvmethod
!! Flag for the matrix-vector method to be used; optional, should be one of
!! OMMP_MATV_ if not provided [[mod_constants:OMMP_MATV_DEFAULT]] is used.
logical, intent(in), optional :: arg_ipd_mask(sys_obj%eel%n_ipd)
!! Logical mask to skip calculation of one of the two set of dipoles
!! in AMOEBA calculations (eg. when MM part's field is not taken into
!! account, both P and D field are just external field, so there is no
!! reason to compute it twice). If n_ipd == 1 this is always considered
!! true.
real(rp), dimension(:, :), allocatable :: e_vec, ipd0
real(rp), dimension(:), allocatable :: inv_diag
integer(ip) :: i, n, solver, mvmethod
logical :: ipd_mask(sys_obj%eel%n_ipd), amoeba
type(ommp_electrostatics_type), pointer :: eel
abstract interface
subroutine mv(eel, x, y, dodiag)
use mod_memory, only: rp, ip
use mod_electrostatics, only : ommp_electrostatics_type
type(ommp_electrostatics_type), intent(in) :: eel
real(rp), dimension(3*eel%pol_atoms), intent(in) :: x
real(rp), dimension(3*eel%pol_atoms), intent(out) :: y
logical, intent(in) :: dodiag
end subroutine mv
end interface
procedure(mv), pointer :: matvec
abstract interface
subroutine pc(eel, x, y)
use mod_memory, only: rp, ip
use mod_electrostatics, only : ommp_electrostatics_type
type(ommp_electrostatics_type), intent(in) :: eel
real(rp), dimension(3*eel%pol_atoms), intent(in) :: x
real(rp), dimension(3*eel%pol_atoms), intent(out) :: y
end subroutine pc
end interface
procedure(pc), pointer :: precond
call time_push()
! Shortcuts
eel => sys_obj%eel
amoeba = sys_obj%amoeba
! Defaults for safety
matvec => TMatVec_incore
precond => PolVec
if(eel%pol_atoms == 0) then
! If the system is not polarizable, there is nothing to do.
call time_pull('Polarization routine')
return
end if
! Handling of optional arguments
if(present(arg_solver)) then
solver = arg_solver
if(solver == OMMP_SOLVER_NONE) solver = eel%def_solver
else
solver = eel%def_solver
end if
if(present(arg_mvmethod)) then
mvmethod = arg_mvmethod
if(mvmethod == OMMP_MATV_NONE) mvmethod = eel%def_matv
else
mvmethod = eel%def_matv
end if
if(present(arg_ipd_mask) .and. eel%n_ipd > 1) then
ipd_mask = arg_ipd_mask
else
ipd_mask = .true.
end if
! Dimension of the system
n = 3*eel%pol_atoms
call mallocate('polarization [ipd0]', n, eel%n_ipd, ipd0)
call mallocate('polarization [e_vec]', n, eel%n_ipd, e_vec)
! Allocate and compute dipole polarization tensor, if needed
if(mvmethod == OMMP_MATV_INCORE .or. &
solver == OMMP_SOLVER_INVERSION) then
if(.not. allocated(eel%tmat)) then !TODO move this in create_tmat
call ommp_message("Allocating T matrix.", OMMP_VERBOSE_DEBUG)
call mallocate('polarization [TMat]',n,n,eel%tmat)
call create_TMat(eel)
end if
end if
! Reshape electric field matrix into a vector
! direct field for Wang and Amoeba
! polarization field just for Amoeba
if(amoeba) then
if(ipd_mask(_amoeba_D_)) &
e_vec(:, _amoeba_D_) = reshape(e(:,:,_amoeba_D_), (/ n /))
if(ipd_mask(_amoeba_P_)) &
e_vec(:, _amoeba_P_) = reshape(e(:,:,_amoeba_P_), (/ n /))
else
e_vec(:, 1) = reshape(e(:,:, 1), (/ n /))
end if
! Initialization of dipoles
ipd0 = 0.0_rp
if(solver /= OMMP_SOLVER_INVERSION) then
! Create a guess for dipoles
if(eel%ipd_use_guess) then
if(amoeba) then
if(ipd_mask(_amoeba_D_)) then
ipd0(:,_amoeba_D_) = &
reshape(eel%ipd(:,:,_amoeba_D_), [n])
end if
if(ipd_mask(_amoeba_P_)) then
ipd0(:, _amoeba_P_) = &
reshape(eel%ipd(:,:,_amoeba_P_), [n])
end if
else
! call PolVec(eel, e_vec(:,1), ipd0(:,1))
ipd0(:, 1) = &
reshape(eel%ipd(:,:,1), [n])
end if
end if
select case(mvmethod)
case(OMMP_MATV_INCORE)
call ommp_message("Matrix-Vector will be performed in-memory", &
OMMP_VERBOSE_HIGH)
matvec => TMatVec_incore
case(OMMP_MATV_DIRECT)
call ommp_message("Matrix-Vector will be performed on-the-fly", &
OMMP_VERBOSE_HIGH)
matvec => TMatVec_otf
case default
call fatal_error("Unknown matrix-vector method requested")
end select
end if
select case (solver)
case(OMMP_SOLVER_CG)
! For now we do not have any other option.
precond => PolVec
if(amoeba) then
if(ipd_mask(_amoeba_D_)) &
call conjugate_gradient_solver(n, &
e_vec(:,_amoeba_D_), &
ipd0(:,_amoeba_D_), &
eel, matvec, precond)
! If both sets have to be computed and there is no input
! guess, just use D as guess for P, not a big gain but still
! something
if(ipd_mask(_amoeba_D_) .and. ipd_mask(_amoeba_P_) &
.and. .not. eel%ipd_use_guess) &
ipd0(:,_amoeba_P_) = ipd0(:,_amoeba_D_)
if(ipd_mask(_amoeba_P_)) &
call conjugate_gradient_solver(n, &
e_vec(:,_amoeba_P_), &
ipd0(:,_amoeba_P_), &
eel, matvec, precond)
else
call conjugate_gradient_solver(n, e_vec(:,1), ipd0(:,1), &
eel, matvec, precond)
end if
case(OMMP_SOLVER_DIIS)
! Create a vector containing inverse of diagonal of T matrix
call mallocate('polarization [inv_diag]', n, inv_diag)
!$omp parallel do default(shared) private(i) schedule(static)
do i=1, eel%pol_atoms
inv_diag(3*(i-1)+1:3*(i-1)+3) = eel%pol(i)
end do
if(amoeba) then
if(ipd_mask(_amoeba_D_)) &
call jacobi_diis_solver(n, &
e_vec(:,_amoeba_D_), &
ipd0(:,_amoeba_D_), &
eel, matvec, inv_diag)
! If both sets have to be computed and there is no input
! guess, just use D as guess for P, not a big gain but still
! something
if(ipd_mask(_amoeba_D_) .and. ipd_mask(_amoeba_P_) &
.and. .not. eel%ipd_use_guess) &
ipd0(:,_amoeba_P_) = ipd0(:,_amoeba_D_)
if(ipd_mask(_amoeba_P_)) &
call jacobi_diis_solver(n, &
e_vec(:,_amoeba_P_), &
ipd0(:,_amoeba_P_), &
eel, matvec, inv_diag)
else
call jacobi_diis_solver(n, e_vec(:,1), ipd0(:,1), &
eel, matvec, inv_diag)
end if
call mfree('polarization [inv_diag]', inv_diag)
case(OMMP_SOLVER_INVERSION)
if(amoeba) then
if(ipd_mask(_amoeba_D_)) &
call inversion_solver(n, &
e_vec(:,_amoeba_D_), &
ipd0(:,_amoeba_D_), eel%TMat)
if(ipd_mask(_amoeba_P_)) &
call inversion_solver(n, &
e_vec(:,_amoeba_P_), &
ipd0(:,_amoeba_P_), eel%TMat)
else
call inversion_solver(n, e_vec(:,1), ipd0(:,1), eel%TMat)
end if
case default
call fatal_error("Unknown solver for calculation of the induced point dipoles")
end select
! Reshape dipole vector into the matrix
eel%ipd = reshape(ipd0, (/3_ip, eel%pol_atoms, eel%n_ipd/))
eel%ipd_done = .true. !! TODO Maybe check convergence...
eel%ipd_use_guess = .true.
call mfree('polarization [ipd0]', ipd0)
call mfree('polarization [e_vec]', e_vec)
call time_pull('Polarization routine')
end subroutine polarization