field_extD2D Subroutine

public subroutine field_extD2D(eel, ext_ipd, E)

Computes the electric field of a trial set of induced point dipoles at polarizable sites. This is intended to be used as matrix-vector routine in the solution of the linear system.

Arguments

Type IntentOptional Attributes Name
type(ommp_electrostatics_type), intent(in) :: eel

Data structure for electrostatic part of the system
real(kind=rp), intent(in) :: ext_ipd(3,eel%pol_atoms)

External induced point dipoles at polarizable sites
real(kind=rp), intent(inout) :: E(3,eel%pol_atoms)

Electric field (results will be added)

Calls

proc~~field_extd2d~~CallsGraph proc~field_extd2d field_extD2D proc~fmm_init fmm_init proc~field_extd2d->proc~fmm_init proc~prepare_fmm_ext_ipd prepare_fmm_ext_ipd proc~field_extd2d->proc~prepare_fmm_ext_ipd proc~cart_propfar_at_ipart cart_propfar_at_ipart proc~field_extd2d->proc~cart_propfar_at_ipart proc~mu_elec_prop mu_elec_prop proc~field_extd2d->proc~mu_elec_prop proc~damped_coulomb_kernel damped_coulomb_kernel proc~field_extd2d->proc~damped_coulomb_kernel proc~ntot_sph_harm ntot_sph_harm proc~fmm_init->proc~ntot_sph_harm proc~prepare_fmmm_constants prepare_fmmm_constants proc~fmm_init->proc~prepare_fmmm_constants interface~mallocate mallocate proc~prepare_fmm_ext_ipd->interface~mallocate proc~fmm_solve_for_multipoles fmm_solve_for_multipoles proc~prepare_fmm_ext_ipd->proc~fmm_solve_for_multipoles interface~mfree mfree proc~prepare_fmm_ext_ipd->interface~mfree proc~cart_propfar_at_ipart->proc~ntot_sph_harm proc~fmm_l2l fmm_l2l proc~cart_propfar_at_ipart->proc~fmm_l2l proc~coulomb_kernel coulomb_kernel proc~damped_coulomb_kernel->proc~coulomb_kernel proc~fatal_error fatal_error proc~damped_coulomb_kernel->proc~fatal_error 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_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~make_vcnk make_vcnk proc~prepare_fmmm_constants->proc~make_vcnk proc~make_m2l_ztranslate_coef make_m2l_ztranslate_coef proc~prepare_fmmm_constants->proc~make_m2l_ztranslate_coef proc~make_vscales make_vscales proc~prepare_fmmm_constants->proc~make_vscales proc~fmm_solve_for_multipoles->proc~fatal_error proc~tree_p2m tree_p2m proc~fmm_solve_for_multipoles->proc~tree_p2m proc~tree_m2m tree_m2m proc~fmm_solve_for_multipoles->proc~tree_m2m proc~tree_l2l tree_l2l proc~fmm_solve_for_multipoles->proc~tree_l2l proc~tree_m2l tree_m2l proc~fmm_solve_for_multipoles->proc~tree_m2l proc~time_pull time_pull proc~fmm_solve_for_multipoles->proc~time_pull proc~time_push time_push proc~fmm_solve_for_multipoles->proc~time_push 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_free3 i_free3 interface~mfree->proc~i_free3 proc~r_free1 r_free1 interface~mfree->proc~r_free1 proc~i_free1 i_free1 interface~mfree->proc~i_free1 proc~fmm_l2l_rotation_work fmm_l2l_rotation_work proc~fmm_l2l->proc~fmm_l2l_rotation_work proc~make_vfact make_vfact proc~fmm_l2l->proc~make_vfact proc~coulomb_kernel->proc~fatal_error proc~ommp_message ommp_message proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~chk_free chk_free proc~i_free2->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~tree_p2m->proc~ntot_sph_harm proc~fmm_error fmm_error proc~tree_p2m->proc~fmm_error proc~fmm_m2m fmm_m2m proc~tree_p2m->proc~fmm_m2m proc~r_free3->proc~chk_free proc~make_m2l_ztranslate_coef->proc~make_vcnk proc~r_free2->proc~chk_free proc~l_free1->proc~chk_free proc~l_free2->proc~chk_free proc~close_output->proc~ommp_message proc~i_free3->proc~chk_free proc~carttosph carttosph proc~fmm_l2l_rotation_work->proc~carttosph proc~fmm_sph_rotate_oxz_work fmm_sph_rotate_oxz_work proc~fmm_l2l_rotation_work->proc~fmm_sph_rotate_oxz_work proc~trgev trgev proc~fmm_l2l_rotation_work->proc~trgev proc~fmm_l2l_ztranslate_work fmm_l2l_ztranslate_work proc~fmm_l2l_rotation_work->proc~fmm_l2l_ztranslate_work proc~fmm_sph_rotate_oz_adj_work fmm_sph_rotate_oz_adj_work proc~fmm_l2l_rotation_work->proc~fmm_sph_rotate_oz_adj_work proc~fmm_sph_rotate_oz_work fmm_sph_rotate_oz_work proc~fmm_l2l_rotation_work->proc~fmm_sph_rotate_oz_work proc~tree_m2m->proc~ntot_sph_harm proc~tree_m2m->proc~fmm_m2m proc~tree_l2l->proc~ntot_sph_harm proc~tree_l2l->proc~fmm_l2l proc~tree_m2l->proc~ntot_sph_harm proc~fmm_m2l fmm_m2l proc~tree_m2l->proc~fmm_m2l proc~time_pull->proc~fatal_error proc~time_pull->proc~ommp_message proc~mem_stat mem_stat proc~time_pull->proc~mem_stat proc~r_free1->proc~chk_free proc~r_alloc1->proc~chk_alloc proc~r_alloc1->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_alloc2->proc~chk_alloc proc~r_alloc2->proc~memory_init proc~l_alloc1->proc~chk_alloc proc~l_alloc1->proc~memory_init proc~i_free1->proc~chk_free proc~time_push->proc~fatal_error proc~time_push->proc~mem_stat proc~mem_stat->proc~memory_init proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error proc~fmm_m2m_rotation_work fmm_m2m_rotation_work proc~fmm_m2m->proc~fmm_m2m_rotation_work proc~fmm_m2l_rotation_work fmm_m2l_rotation_work proc~fmm_m2l->proc~fmm_m2l_rotation_work proc~fmm_m2m_rotation_work->proc~carttosph proc~fmm_m2m_rotation_work->proc~fmm_sph_rotate_oxz_work proc~fmm_m2m_rotation_work->proc~trgev proc~fmm_m2m_rotation_work->proc~fmm_sph_rotate_oz_adj_work proc~fmm_m2m_rotation_work->proc~fmm_sph_rotate_oz_work proc~fmm_m2m_ztranslate_work fmm_m2m_ztranslate_work proc~fmm_m2m_rotation_work->proc~fmm_m2m_ztranslate_work proc~fmm_m2l_rotation_work->proc~carttosph proc~fmm_m2l_rotation_work->proc~fmm_sph_rotate_oxz_work proc~fmm_m2l_rotation_work->proc~trgev proc~fmm_m2l_rotation_work->proc~fmm_sph_rotate_oz_adj_work proc~fmm_m2l_rotation_work->proc~fmm_sph_rotate_oz_work proc~fmm_m2l_ztranslate_work fmm_m2l_ztranslate_work proc~fmm_m2l_rotation_work->proc~fmm_m2l_ztranslate_work

Called by

proc~~field_extd2d~~CalledByGraph proc~field_extd2d field_extD2D proc~tmatvec_otf TMatVec_otf proc~tmatvec_otf->proc~field_extd2d

Contents

Source Code

field_extD2D

Source Code

    subroutine field_extD2D(eel, ext_ipd, E)
        !! Computes the electric field of a trial set of induced point dipoles
        !! at polarizable sites. This is intended to be used as matrix-vector
        !! routine in the solution of the linear system.
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Data structure for electrostatic part of the system
        real(rp), intent(in) :: ext_ipd(3, eel%pol_atoms)
        !! External induced point dipoles at polarizable sites
        real(rp), intent(inout) :: E(3, eel%pol_atoms)
        !! Electric field (results will be added)

        integer(ip) :: i, j, ipol, jpol, ij, idx
        logical :: to_scale, to_do
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), scalf
        type(fmm_type), allocatable :: fmm_ipd

        if(eel%use_fmm) then
            allocate(fmm_ipd)
            call fmm_init(fmm_ipd, eel%fmm_maxl_pol, eel%tree)
            call prepare_fmm_ext_ipd(eel, fmm_ipd, ext_ipd)

            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,ij,ipol,jpol,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
            do ipol=1, eel%pol_atoms 
                i = eel%polar_mm(ipol)
                tmpE = 0.0
                call cart_propfar_at_ipart(fmm_ipd, i, &
                                           .false., tmpV, &
                                           .true. , tmpE, &
                                           .false., tmpEgr, &
                                           .false., tmpHE)
                E(:, ipol) = tmpE
            end do
            
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,ij,ipol,jpol,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
            do ipol=1, eel%pol_atoms 
                i = eel%polar_mm(ipol)
                 
                ! Near field is computed internally because dumped kernel is required
                do ij=eel%fmm_near_field_list%ri(i), eel%fmm_near_field_list%ri(i+1)-1
                    j = eel%fmm_near_field_list%ci(ij)
                    jpol = eel%mm_polar(j)
                    ! If the atom is not polarizable, skip
                    if(jpol < 1) cycle 

                    !loop on target
                    to_do = .true.
                    to_scale = .false.
                    scalf = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_P_P%ri(ipol), eel%list_P_P%ri(ipol+1)-1
                        if(eel%list_P_P%ci(idx) == jpol) then
                            to_scale = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale) then
                        to_do = eel%todo_P_P(idx)
                        scalf = eel%scalef_P_P(idx)
                    end if
                    
                    if(to_do) then
                        call damped_coulomb_kernel(eel, j, i,& 
                                                   2, kernel(1:3), dr)
                        
                        tmpE = 0.0_rp

                        call mu_elec_prop(ext_ipd(:,jpol), dr, kernel, .false., tmpV, &
                                        .true., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)
                        if(to_scale) then
                            E(:, ipol) = E(:, ipol) + tmpE * scalf
                        else
                            E(:, ipol) = E(:, ipol) + tmpE
                        end if
                    end if
                end do
            end do
            
            if(allocated(eel%list_P_P_fmm_far)) then
                ! Now remove screened interactions from far-field
                do ipol=1, eel%pol_atoms
                    i = eel%polar_mm(ipol)
                    
                    do idx=eel%list_P_P_fmm_far%ri(ipol), eel%list_P_P_fmm_far%ri(ipol+1)-1
                        jpol = eel%list_P_P_fmm_far%ci(idx)
                        j = eel%polar_mm(jpol)

                        scalf = 1.0 - eel%scalef_P_P_fmm_far(idx)
                        
                        call damped_coulomb_kernel(eel, j, i,& 
                                                    2, kernel(1:3), dr)
                        
                        tmpE = 0.0_rp
                        call mu_elec_prop(ext_ipd(:,jpol), dr, kernel, .false., tmpV, &
                                        .true., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)
                        
                        E(:, ipol) = E(:, ipol) - tmpE * scalf
                    end do
                end do
            end if
            deallocate(fmm_ipd%multipoles)
            deallocate(fmm_ipd%local_expansion)
            deallocate(fmm_ipd)
        else
        
        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,to_do,to_scale,scalf,idx,tmpV,tmpE,tmpEgr,tmpHE,kernel,dr)
        do j=1, eel%pol_atoms
            do i=1, eel%pol_atoms
                if(j == i) cycle
                !loop on target
                to_do = .true.
                to_scale = .false.
                scalf = 1.0

                ! Check if the element should be scaled
                do idx=eel%list_P_P%ri(i), eel%list_P_P%ri(i+1)-1
                    if(eel%list_P_P%ci(idx) == j) then
                        to_scale = .true.
                        exit
                    end if
                end do

                !If it should set the correct variables
                if(to_scale) then
                    to_do = eel%todo_P_P(idx)
                    scalf = eel%scalef_P_P(idx)
                end if
                
                if(to_do) then
                    call damped_coulomb_kernel(eel, eel%polar_mm(i), &
                                               eel%polar_mm(j),& 
                                               2, kernel(1:3), dr)
                    
                    tmpE = 0.0_rp

                    call mu_elec_prop(ext_ipd(:,i), dr, kernel, .false., tmpV, &
                                      .true., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)
                    if(to_scale) then
                        E(:, j) = E(:, j) + tmpE * scalf
                    else
                        E(:, j) = E(:, j) + tmpE
                    end if
                end if
            end do
        end do
        end if
    end subroutine field_extD2D