ddx_cosmo.f90

Go to the documentation of this file.

 
 
module ddx_cosmo
! Get ddx-operators
use ddx_operators
use ddx_multipolar_solutes
implicit none
 
 
contains
 
subroutine cosmo_energy(constants, state, esolv, ddx_error)
    implicit none
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_state_type), intent(in) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
    real(dp), intent(out) :: esolv
    real(dp), external :: ddot
 
    ! dummy operation on unused interface arguments
    if (ddx_error % flag .eq. 0) continue
 
    esolv = pt5*ddot(constants % n, state % xs, 1, state % psi, 1)
end subroutine cosmo_energy
 
subroutine cosmo_setup(params, constants, workspace, state, phi_cav, psi, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
    real(dp), intent(in) :: phi_cav(constants % ncav)
    real(dp), intent(in) :: psi(constants % nbasis, params % nsph)
 
    ! dummy operation on unused interface arguments
    if (ddx_error % flag .eq. 0) continue
 
    call cav_to_spherical(params, constants, workspace, phi_cav, &
        & state % phi)
    state % phi = - state % phi
    state % phi_cav = phi_cav
    state % psi = psi
    state % rhs_done = .true.
    state % adjoint_rhs_done = .true.
end subroutine cosmo_setup
 
subroutine cosmo_guess(params, constants, workspace, state, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
 
    ! apply the diagonal preconditioner as a guess
    call ldm1x(params, constants, workspace, state % phi, state % xs, ddx_error)
 
end subroutine cosmo_guess
 
subroutine cosmo_guess_adjoint(params, constants, workspace, state, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
 
    ! apply the diagonal preconditioner as a guess
    call ldm1x(params, constants, workspace, state % psi, state % s, ddx_error)
 
end subroutine cosmo_guess_adjoint
 
subroutine cosmo_solve(params, constants, workspace, state, tol, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    real(dp), intent(in) :: tol
    type(ddx_error_type), intent(inout) :: ddx_error
    ! local variables
    real(dp) :: start_time, finish_time
 
    state % xs_niter =  params % maxiter
    start_time = omp_get_wtime()
    call jacobi_diis(params, constants, workspace, tol, state % phi, &
        & state % xs, state % xs_niter, state % xs_rel_diff, lx, ldm1x, &
        & hnorm, ddx_error)
    finish_time = omp_get_wtime()
    state % xs_time = finish_time - start_time
 
end subroutine cosmo_solve
 
subroutine cosmo_solve_adjoint(params, constants, workspace, state, tol, &
        & ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    real(dp), intent(in) :: tol
    type(ddx_error_type), intent(inout) :: ddx_error
    ! local variables
    real(dp) :: start_time, finish_time
 
    state % s_niter = params % maxiter
    start_time = omp_get_wtime()
    call jacobi_diis(params, constants, workspace, tol, state % psi, &
        & state % s, state % s_niter, state % s_rel_diff, lstarx, ldm1x, &
        & hnorm, ddx_error)
    finish_time = omp_get_wtime()
    state % s_time = finish_time - start_time
 
    state % q = state % s
 
    call cosmo_derivative_setup(params, constants, workspace, state)
 
end subroutine cosmo_solve_adjoint
 
subroutine cosmo_solvation_force_terms(params, constants, workspace, &
        & state, e_cav, force, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
    real(dp), intent(in) :: e_cav(3, constants % ncav)
    real(dp), intent(inout) :: force(3, params % nsph)
    ! local variables
    real(dp) :: start_time, finish_time
    integer :: isph
 
    ! dummy operation on unused interface arguments
    if (ddx_error % flag .eq. 0) continue
 
    start_time = omp_get_wtime()
 
    force = zero
    do isph = 1, params % nsph
        call contract_grad_l(params, constants, isph, state % xs, &
            & state % sgrid, workspace % tmp_vylm(:, 1), &
            & workspace % tmp_vdylm(:, :, 1), workspace % tmp_vplm(:, 1), &
            & workspace % tmp_vcos(:, 1), workspace % tmp_vsin(:, 1), &
            & force(:, isph))
        call contract_grad_u(params, constants, isph, state % sgrid, &
            & state % phi_grid, force(:, isph))
    end do
 
    force = -pt5*force
 
    call zeta_grad(params, constants, state, e_cav, force)
 
    finish_time = omp_get_wtime()
    state % force_time = finish_time - start_time
 
end subroutine cosmo_solvation_force_terms
 
subroutine cosmo_solvation_force_dr_terms(params, constants, workspace, &
        & state, e_cav, force, dr, ddx_error)
    implicit none
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
    type(ddx_error_type), intent(inout) :: ddx_error
    real(dp), intent(in) :: e_cav(3, constants % ncav)
    real(dp), intent(inout) :: force(3, params % nsph)
    real(dp), intent(inout) :: dr(params % nsph)
    ! local variables
    real(dp) :: start_time, finish_time
    integer :: isph
 
    ! dummy operation on unused interface arguments
    if (ddx_error % flag .eq. 0) continue
 
    start_time = omp_get_wtime()
 
    force = zero
    dr = zero
    do isph = 1, params % nsph
        call contract_grad_l(params, constants, isph, state % xs, &
            & state % sgrid, workspace % tmp_vylm(:, 1), &
            & workspace % tmp_vdylm(:, :, 1), workspace % tmp_vplm(:, 1), &
            & workspace % tmp_vcos(:, 1), workspace % tmp_vsin(:, 1), &
            & force(:, isph), dr(isph))
        call contract_grad_u(params, constants, isph, state % sgrid, &
            & state % phi_grid, force(:, isph), dr(isph))
    end do
 
    force = -pt5*force
    dr = -pt5*dr
 
    call zeta_grad(params, constants, state, e_cav, force)
    call zeta_grad_dr(params, constants, state, e_cav, dr)
 
    finish_time = omp_get_wtime()
    state % force_time = finish_time - start_time
 
end subroutine cosmo_solvation_force_dr_terms
 
 
subroutine cosmo_derivative_setup(params, constants, workspace, state)
    type(ddx_params_type), intent(in) :: params
    type(ddx_constants_type), intent(in) :: constants
    type(ddx_workspace_type), intent(inout) :: workspace
    type(ddx_state_type), intent(inout) :: state
 
    real(dp), external :: ddot
    integer :: icav, isph, igrid
 
    ! dummy operation on unused interface arguments
    if (allocated(workspace % tmp_pot)) continue
 
    ! Get values of S on the grid
    call ddeval_grid_work(constants % nbasis, params % ngrid, params % nsph, &
        & constants % vgrid, constants % vgrid_nbasis, one, state % s, zero, &
        & state % sgrid)
    ! Get the values of phi on the grid
    call ddcav_to_grid_work(params % ngrid, params % nsph, constants % ncav, &
        & constants % icav_ia, constants % icav_ja, state % phi_cav, &
        & state % phi_grid)
 
    ! assemble the intermediate zeta: S weighted by U evaluated on the
    ! exposed grid points.
    icav = 0
    do isph = 1, params % nsph
        do igrid = 1, params % ngrid
            if (constants % ui(igrid, isph) .ne. zero) then
                icav = icav + 1
                state % zeta(icav) = constants % wgrid(igrid) * &
                    & constants % ui(igrid, isph) * ddot(constants % nbasis, &
                    & constants % vgrid(1, igrid), 1, state % s(1, isph), 1)
            end if
        end do
    end do
 
end subroutine cosmo_derivative_setup
 
end module ddx_cosmo