polarization Subroutine

public subroutine polarization(sys_obj, e, arg_solver, arg_mvmethod, arg_ipd_mask)

Uses

  • proc~~polarization~~UsesGraph proc~polarization polarization module~mod_profiling mod_profiling proc~polarization->module~mod_profiling module~mod_constants mod_constants proc~polarization->module~mod_constants module~mod_io mod_io proc~polarization->module~mod_io module~mod_memory mod_memory proc~polarization->module~mod_memory module~mod_solvers mod_solvers proc~polarization->module~mod_solvers module~mod_profiling->module~mod_constants module~mod_profiling->module~mod_io module~mod_profiling->module~mod_memory iso_c_binding iso_c_binding module~mod_constants->iso_c_binding module~mod_io->module~mod_constants module~mod_memory->module~mod_constants module~mod_memory->module~mod_io module~mod_memory->iso_c_binding module~mod_solvers->module~mod_constants module~mod_solvers->module~mod_io module~mod_solvers->module~mod_memory module~mod_electrostatics mod_electrostatics module~mod_solvers->module~mod_electrostatics module~mod_electrostatics->module~mod_profiling module~mod_electrostatics->module~mod_constants module~mod_electrostatics->module~mod_io module~mod_electrostatics->module~mod_memory 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~fmmlib_interface fmmlib_interface module~mod_electrostatics->module~fmmlib_interface module~mod_adjacency_mat->module~mod_memory module~mod_topology->module~mod_memory module~mod_topology->module~mod_adjacency_mat module~fmmlib_interface->module~mod_constants module~mod_tree mod_tree module~fmmlib_interface->module~mod_tree module~mod_ribtree mod_ribtree module~fmmlib_interface->module~mod_ribtree module~mod_harmonics mod_harmonics module~fmmlib_interface->module~mod_harmonics module~mod_fmm_utils mod_fmm_utils module~fmmlib_interface->module~mod_fmm_utils module~mod_fmm mod_fmm module~fmmlib_interface->module~mod_fmm module~mod_octatree mod_octatree module~fmmlib_interface->module~mod_octatree module~mod_tree->module~mod_constants module~mod_tree->module~mod_adjacency_mat module~mod_tree->module~mod_fmm_utils module~mod_ribtree->module~mod_profiling module~mod_ribtree->module~mod_constants module~mod_ribtree->module~mod_tree module~mod_ribtree->module~mod_fmm_utils module~mod_harmonics->module~mod_constants module~mod_harmonics->module~mod_fmm_utils module~mod_fmm_utils->module~mod_constants module~mod_fmm->module~mod_constants module~mod_fmm->module~mod_tree module~mod_fmm->module~mod_harmonics module~mod_fmm->module~mod_fmm_utils module~mod_octatree->module~mod_profiling module~mod_octatree->module~mod_constants module~mod_octatree->module~mod_tree module~mod_octatree->module~mod_fmm_utils

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...

Arguments

Type IntentOptional 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.

Calls

proc~~polarization~~CallsGraph proc~polarization polarization proc~time_push time_push proc~polarization->proc~time_push interface~mallocate mallocate proc~polarization->interface~mallocate proc~create_tmat create_tmat proc~polarization->proc~create_tmat proc~time_pull time_pull proc~polarization->proc~time_pull proc~ommp_message ommp_message proc~polarization->proc~ommp_message proc~fatal_error fatal_error proc~polarization->proc~fatal_error proc~conjugate_gradient_solver conjugate_gradient_solver proc~polarization->proc~conjugate_gradient_solver proc~jacobi_diis_solver jacobi_diis_solver proc~polarization->proc~jacobi_diis_solver interface~mfree mfree proc~polarization->interface~mfree proc~inversion_solver inversion_solver proc~polarization->proc~inversion_solver proc~time_push->proc~fatal_error proc~mem_stat mem_stat proc~time_push->proc~mem_stat proc~i_alloc1 i_alloc1 interface~mallocate->proc~i_alloc1 proc~i_alloc3 i_alloc3 interface~mallocate->proc~i_alloc3 proc~l_alloc2 l_alloc2 interface~mallocate->proc~l_alloc2 proc~r_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~r_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~create_tmat->proc~ommp_message proc~dipole_t dipole_T proc~create_tmat->proc~dipole_t proc~time_pull->proc~ommp_message proc~time_pull->proc~fatal_error proc~time_pull->proc~mem_stat proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~conjugate_gradient_solver->interface~mallocate proc~conjugate_gradient_solver->proc~ommp_message proc~conjugate_gradient_solver->proc~fatal_error proc~conjugate_gradient_solver->interface~mfree proc~jacobi_diis_solver->interface~mallocate proc~jacobi_diis_solver->proc~ommp_message proc~jacobi_diis_solver->proc~fatal_error proc~jacobi_diis_solver->interface~mfree proc~diis diis proc~jacobi_diis_solver->proc~diis proc~rmsvec rmsvec proc~jacobi_diis_solver->proc~rmsvec proc~i_free2 i_free2 interface~mfree->proc~i_free2 proc~r_free3 r_free3 interface~mfree->proc~r_free3 proc~r_free2 r_free2 interface~mfree->proc~r_free2 proc~l_free1 l_free1 interface~mfree->proc~l_free1 proc~l_free2 l_free2 interface~mfree->proc~l_free2 proc~i_free1 i_free1 interface~mfree->proc~i_free1 proc~i_free3 i_free3 interface~mfree->proc~i_free3 proc~r_free1 r_free1 interface~mfree->proc~r_free1 proc~inversion_solver->interface~mallocate proc~inversion_solver->interface~mfree dgemm dgemm proc~inversion_solver->dgemm dgetri dgetri proc~inversion_solver->dgetri dgetrf dgetrf proc~inversion_solver->dgetrf proc~chk_free chk_free proc~i_free2->proc~chk_free proc~r_free3->proc~chk_free proc~chk_alloc chk_alloc proc~i_alloc1->proc~chk_alloc proc~memory_init memory_init proc~i_alloc1->proc~memory_init proc~i_alloc3->proc~chk_alloc proc~i_alloc3->proc~memory_init proc~l_alloc2->proc~chk_alloc proc~l_alloc2->proc~memory_init proc~screening_rules screening_rules proc~dipole_t->proc~screening_rules proc~damped_coulomb_kernel damped_coulomb_kernel proc~dipole_t->proc~damped_coulomb_kernel proc~r_free2->proc~chk_free proc~close_output->proc~ommp_message proc~l_free1->proc~chk_free proc~l_free2->proc~chk_free proc~diis->interface~mallocate proc~diis->interface~mfree dgesv dgesv proc~diis->dgesv proc~makeb makeb proc~diis->proc~makeb proc~i_free1->proc~chk_free proc~i_free3->proc~chk_free proc~r_free1->proc~chk_free proc~mem_stat->proc~memory_init proc~r_alloc3->proc~chk_alloc proc~r_alloc3->proc~memory_init proc~i_alloc2->proc~chk_alloc proc~i_alloc2->proc~memory_init proc~r_alloc1->proc~chk_alloc proc~r_alloc1->proc~memory_init proc~r_alloc2->proc~chk_alloc proc~r_alloc2->proc~memory_init proc~l_alloc1->proc~chk_alloc proc~l_alloc1->proc~memory_init proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error proc~screening_rules->proc~fatal_error proc~damped_coulomb_kernel->proc~fatal_error proc~coulomb_kernel coulomb_kernel proc~damped_coulomb_kernel->proc~coulomb_kernel proc~coulomb_kernel->proc~fatal_error

Called by

proc~~polarization~~CalledByGraph proc~polarization polarization proc~polelec_geomgrad polelec_geomgrad proc~polelec_geomgrad->proc~polarization proc~ommp_get_polelec_energy ommp_get_polelec_energy proc~ommp_get_polelec_energy->proc~polarization proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~polarization proc~ommp_polelec_geomgrad ommp_polelec_geomgrad proc~ommp_polelec_geomgrad->proc~polelec_geomgrad proc~c_ommp_get_polelec_energy C_ommp_get_polelec_energy proc~c_ommp_get_polelec_energy->proc~ommp_get_polelec_energy proc~ommp_get_full_ele_energy ommp_get_full_ele_energy proc~ommp_get_full_ele_energy->proc~ommp_get_polelec_energy proc~ommp_set_external_field_nomm ommp_set_external_field_nomm proc~ommp_set_external_field_nomm->proc~ommp_set_external_field proc~ommp_full_geomgrad ommp_full_geomgrad proc~ommp_full_geomgrad->proc~polelec_geomgrad proc~c_ommp_set_external_field_nomm C_ommp_set_external_field_nomm proc~c_ommp_set_external_field_nomm->proc~ommp_set_external_field proc~c_ommp_set_external_field C_ommp_set_external_field proc~c_ommp_set_external_field->proc~ommp_set_external_field proc~c_ommp_polelec_geomgrad C_ommp_polelec_geomgrad proc~c_ommp_polelec_geomgrad->proc~ommp_polelec_geomgrad proc~c_ommp_full_geomgrad C_ommp_full_geomgrad proc~c_ommp_full_geomgrad->proc~ommp_full_geomgrad proc~c_ommp_get_full_ele_energy C_ommp_get_full_ele_energy proc~c_ommp_get_full_ele_energy->proc~ommp_get_full_ele_energy proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_ele_energy proc~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code

polarization

Source Code

    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