damped_coulomb_kernel Subroutine

public subroutine damped_coulomb_kernel(eel, i, j, maxder, res, dr)

Uses

  • proc~~damped_coulomb_kernel~~UsesGraph proc~damped_coulomb_kernel damped_coulomb_kernel module~mod_constants mod_constants proc~damped_coulomb_kernel->module~mod_constants module~mod_memory mod_memory proc~damped_coulomb_kernel->module~mod_memory iso_c_binding iso_c_binding module~mod_constants->iso_c_binding module~mod_memory->module~mod_constants module~mod_memory->iso_c_binding module~mod_io mod_io module~mod_memory->module~mod_io module~mod_io->module~mod_constants

This subroutine computes the damped coulomb kernel between two atoms. Note that this only makes sense between two MM atoms, as it is only used to compute the field that induces the point dipoles!

Arguments

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

Electrostatics data structure
integer(kind=ip), intent(in) :: i

Index of atoms (in MM atom list) for which the kernel should be computed
integer(kind=ip), intent(in) :: j

Index of atoms (in MM atom list) for which the kernel should be computed
integer(kind=ip), intent(in) :: maxder

Maximum derivative to be computed
real(kind=rp), intent(out), dimension(maxder+1) :: res

Results vector
real(kind=rp), intent(out), dimension(3) :: dr

Distance vector between i and j

Calls

proc~~damped_coulomb_kernel~~CallsGraph proc~damped_coulomb_kernel damped_coulomb_kernel 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~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~close_output->proc~ommp_message

Called by

proc~~damped_coulomb_kernel~~CalledByGraph proc~damped_coulomb_kernel damped_coulomb_kernel proc~dipole_t dipole_T proc~dipole_t->proc~damped_coulomb_kernel proc~elec_prop_d2m elec_prop_D2M proc~elec_prop_d2m->proc~damped_coulomb_kernel proc~field_extd2d field_extD2D proc~field_extd2d->proc~damped_coulomb_kernel proc~elec_prop_m2d elec_prop_M2D proc~elec_prop_m2d->proc~damped_coulomb_kernel proc~elec_prop_d2d elec_prop_D2D proc~elec_prop_d2d->proc~damped_coulomb_kernel proc~create_tmat create_tmat proc~create_tmat->proc~dipole_t proc~prepare_polelec prepare_polelec proc~prepare_polelec->proc~elec_prop_d2m proc~prepare_polelec->proc~elec_prop_m2d proc~prepare_polelec->proc~elec_prop_d2d proc~tmatvec_otf TMatVec_otf proc~tmatvec_otf->proc~field_extd2d proc~polarization polarization proc~polarization->proc~create_tmat proc~ommp_get_polelec_energy ommp_get_polelec_energy proc~ommp_get_polelec_energy->proc~prepare_polelec proc~ommp_get_polelec_energy->proc~polarization proc~energy_mm_pol energy_MM_pol proc~ommp_get_polelec_energy->proc~energy_mm_pol proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~prepare_polelec proc~ommp_set_external_field->proc~polarization proc~polelec_geomgrad polelec_geomgrad proc~polelec_geomgrad->proc~prepare_polelec proc~polelec_geomgrad->proc~polarization proc~energy_mm_pol->proc~prepare_polelec proc~ommp_get_full_ele_energy ommp_get_full_ele_energy proc~ommp_get_full_ele_energy->proc~ommp_get_polelec_energy proc~c_ommp_get_polelec_energy C_ommp_get_polelec_energy proc~c_ommp_get_polelec_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~c_ommp_set_external_field C_ommp_set_external_field proc~c_ommp_set_external_field->proc~ommp_set_external_field 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~ommp_polelec_geomgrad ommp_polelec_geomgrad proc~ommp_polelec_geomgrad->proc~polelec_geomgrad proc~ommp_full_geomgrad ommp_full_geomgrad proc~ommp_full_geomgrad->proc~polelec_geomgrad proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_ele_energy 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~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code

damped_coulomb_kernel

Source Code

    subroutine damped_coulomb_kernel(eel, i, j, maxder, res, dr)
        !! This subroutine computes the damped coulomb kernel between two atoms.
        !! Note that this only makes sense between two MM atoms, as it is only used
        !! to compute the field that induces the point dipoles!

        use mod_memory, only: ip, rp
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        integer(ip), intent(in) :: i, j
        !! Index of atoms (in MM atom list) for which the kernel should be
        !! computed
        integer(ip), intent(in) :: maxder
        !! Maximum derivative to be computed
        real(rp), intent(out), dimension(maxder+1) :: res
        !! Results vector
        real(rp), intent(out), dimension(3) :: dr
        !! Distance vector between i and j
        
        real(rp) :: s, u, u3, u4, fexp, eexp
        
        s = eel%thole(i) * eel%thole(j)
        
        ! Compute undamped kernels
        dr = eel%top%cmm(:,j) - eel%top%cmm(:,i)
        call coulomb_kernel(dr, maxder, res)

        if(abs(s) < eps_rp) then
            ! either thole(i) or thole(j) are zero, so the damped kernel
            ! is just the same as the regular one. Job is done.
            return
        end if

        u = 1.0_rp / (res(1) * s)
        u3 = u**3

        if(eel%amoeba .and. u3 * 0.39_rp < 50.0_rp) then
            ! 1.0 / (res(1) * s)**3 * 0.39_rp < 50_rp TODO Remove 0.39, this is a
            ! FF constants, ask someone why this condition is here.
            ! Here basically we multiply the standard interactions kernel for
            ! dumping coefficients. The equations implemented here correspond to 
            ! eq. (5) of 10.1021/jp027815
            fexp = -0.39_rp * u3
            eexp = exp(fexp)
            if(maxder >= 1) res(2) = res(2) * (1.0_rp - eexp)
            if(maxder >= 2) res(3) = res(3) * (1.0_rp - (1.0_rp - fexp) * eexp)
            if(maxder >= 3) res(4) = res(4) * &
                            (1.0_rp - (1.0_rp - fexp + 0.6_rp * fexp * fexp) * eexp)
            if(maxder >= 4) res(5) = res(5) * &
                            (1.0_rp - (1.0_rp - fexp + 18.0_rp/35.0_rp * fexp *  fexp - &
                            9.0_rp/35.0_rp * fexp * fexp * fexp) * eexp)
            if(maxder >= 5) then
                call fatal_error("Damped Coulomb kernel (AMOEBA) only supports up to the 5th derivative")
            end if
        else if(.not. eel%amoeba .and. res(1) > 1_rp/s) then
            ! TODO Again it is not clear to me why condition res(1) > 1_rp/s is here.
            u4 = u3*u
            if(maxder >= 1) res(2) = res(2) * (4.0_rp * u3 - 3.0_rp * u4)
            if(maxder >= 2) res(3) = res(3) * u4
            if(maxder >= 3) res(4) = res(4) * 0.2_rp * u4
            if(maxder >= 4) then
                call fatal_error("Damped Coulomb kernel (WANG) only supports up to the 4th derivative")
            end if
        end if
    end subroutine damped_coulomb_kernel