screening_rules Function

public function screening_rules(eel, i, kind_i, j, kind_j, in_field) result(scalf)

Uses

  • proc~~screening_rules~~UsesGraph proc~screening_rules screening_rules module~mod_constants mod_constants proc~screening_rules->module~mod_constants iso_c_binding iso_c_binding module~mod_constants->iso_c_binding

Utility function used to decide if an interaction between sites i and j should be computed and eventually scaled by a factor. This function is intended to be used in code, for linear scaling code lists should be built. This is written to minimize code repetitions, all the screening rules are handled in two possible cases: 1. rules based on adjacency matrix 2. rules based on AMOEBA polarization groups

Arguments

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

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

Index of source site (MM index is used for static sites, while Pol index is used for polarizable sites)
character, intent(in) :: kind_i

Type of sites i and j in the interaction for which the screening rules are required; possible choices are 'P' (polarizable site) or 'S' (static site). Any other option will cause a fatal error.
integer(kind=ip), intent(in) :: j

Index of target site (MM index is used for static sites, while Pol index is used for polarizable sites)
character, intent(in) :: kind_j

Type of sites i and j in the interaction for which the screening rules are required; possible choices are 'P' (polarizable site) or 'S' (static site). Any other option will cause a fatal error.
character, intent(in) :: in_field

Which screening rules have to be applied? 'D' = screening rules for direct field; 'P' = screening rules for polarization field

Return Value real(kind=rp)

Scale factor for the interaction


Calls

proc~~screening_rules~~CallsGraph proc~screening_rules screening_rules proc~fatal_error fatal_error proc~screening_rules->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~~screening_rules~~CalledByGraph proc~screening_rules screening_rules proc~make_screening_lists make_screening_lists proc~make_screening_lists->proc~screening_rules proc~dipole_t dipole_T proc~dipole_t->proc~screening_rules proc~mmpol_prepare mmpol_prepare proc~mmpol_prepare->proc~make_screening_lists proc~create_tmat create_tmat proc~create_tmat->proc~dipole_t proc~mmpol_init_from_mmp mmpol_init_from_mmp proc~mmpol_init_from_mmp->proc~mmpol_prepare proc~polarization polarization proc~polarization->proc~create_tmat proc~mmpol_init_from_xyz mmpol_init_from_xyz proc~mmpol_init_from_xyz->proc~mmpol_prepare proc~ommp_system_from_qm_helper ommp_system_from_qm_helper proc~ommp_system_from_qm_helper->proc~mmpol_prepare proc~ommp_init_mmp ommp_init_mmp proc~ommp_init_mmp->proc~mmpol_init_from_mmp proc~ommp_init_xyz ommp_init_xyz proc~ommp_init_xyz->proc~mmpol_init_from_xyz proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~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~c_ommp_system_from_qm_helper C_ommp_system_from_qm_helper proc~c_ommp_system_from_qm_helper->proc~ommp_system_from_qm_helper program~test_si_geomgrad test_SI_geomgrad program~test_si_geomgrad->proc~ommp_system_from_qm_helper program~test_si_geomgrad_num test_SI_geomgrad_num program~test_si_geomgrad_num->proc~ommp_system_from_qm_helper proc~c_ommp_init_mmp C_ommp_init_mmp proc~c_ommp_init_mmp->proc~ommp_init_mmp program~test_si_potential test_SI_potential program~test_si_potential->proc~ommp_set_external_field program~test_si_potential->proc~ommp_get_polelec_energy 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_init_xyz C_ommp_init_xyz proc~c_ommp_init_xyz->proc~ommp_init_xyz 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_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_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~c_ommp_polelec_geomgrad C_ommp_polelec_geomgrad proc~c_ommp_polelec_geomgrad->proc~ommp_polelec_geomgrad proc~ommptest_totalqmmm_geomgrad ommptest_totalqmmm_geomgrad proc~ommptest_totalqmmm_geomgrad->proc~ommp_full_geomgrad proc~ommptest_fakeqm_internal_geomgrad ommptest_fakeqm_internal_geomgrad proc~ommptest_fakeqm_internal_geomgrad->proc~ommp_full_geomgrad proc~c_ommp_full_geomgrad C_ommp_full_geomgrad proc~c_ommp_full_geomgrad->proc~ommp_full_geomgrad proc~ommptest_fakeqm_linkatom_geomgrad ommptest_fakeqm_linkatom_geomgrad proc~ommptest_fakeqm_linkatom_geomgrad->proc~ommp_full_geomgrad 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_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

screening_rules

Source Code

    function screening_rules(eel, i, kind_i, j, kind_j, in_field) result(scalf)
        !! Utility function used to decide if an interaction between sites i and j
        !! should be computed and eventually scaled by a factor.
        !! This function is intended to be used in \(\mathcalO(n^2)\) code, for
        !! linear scaling code lists should be built.
        !! This is written to minimize code repetitions, all the screening rules
        !! are handled in two possible cases: 
        !! 1. rules based on adjacency matrix
        !! 2. rules based on AMOEBA polarization groups

        use mod_constants, only: eps_rp

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics object
        integer(ip), intent(in) :: i
        !! Index of source site (MM index is used for static sites, while
        !! Pol index is used for polarizable sites)
        integer(ip), intent(in) :: j
        !! Index of target site (MM index is used for static sites, while
        !! Pol index is used for polarizable sites)
        character, intent(in) :: in_field
        !! Which screening rules have to be applied? 'D' = screening rules
        !! for direct field; 'P' = screening rules for polarization field
        character, intent(in) :: kind_i, kind_j
        !! Type of sites i and j in the interaction for which the screening
        !! rules are required; possible choices are 'P' (polarizable site) or 
        !! 'S' (static site). Any other option will cause a fatal error.
        real(rp) :: scalf
        !! Scale factor for the interaction

        integer(ip) :: ineigh, grp, igrp, pg_i
        integer(ip) :: j_mm, i_mm
        character :: field, interaction
        logical :: amoeba
        type(ommp_topology_type), pointer :: top 
        real(rp) :: myscale(4), myscale_intra(4)

        ! Some shortcut
        top => eel%top
        amoeba = eel%amoeba

        ! Decide which kind of rule should be used
        if(kind_i == 'P') then
            i_mm = eel%polar_mm(i)
        else if(kind_i == 'S') then
            i_mm = i
        else
            call fatal_error('Unexpected value of kind_i in screening_rules')
            i_mm = 0
        end if

        if(kind_j == 'P') then
            j_mm = eel%polar_mm(j)
        else if(kind_j == 'S') then
            j_mm = j
        else
            call fatal_error('Unexpected value of kind_j in screening_rules')
            j_mm = 0
        end if

        myscale = 0.0_rp
        if(kind_j == 'P' .and. kind_i == 'P') then
            ! Use IPD-IPD screening rules
            myscale = eel%uscale
            interaction = 'P' !Pol-pol interaction is named P
        else if(kind_j == 'S' .and. kind_i == 'S') then
            ! Use static multipoles-static multipoles screening rules
            myscale = eel%mscale
            field = '-'
            interaction = 'S' !Static-static interaction is named S
        else
            ! Use static multipoles-IPD screening rules
            if(in_field == 'P' .and. amoeba) then
                myscale = eel%pscale
                myscale_intra = eel%pscale_intra
                field = 'P'
            else if(in_field == 'D' .and. amoeba) then
                myscale = eel%dscale
                field = 'D'
            else if(.not. amoeba) then
                myscale = eel%pscale
                field = 'P'
            else
                call fatal_error('Unexpected value of field in screening rules')
            end if
            interaction = 'M' !Mixed interaction is named M
        end if

        ! Default return value
        scalf = 1.0_rp
        
        if((.not. amoeba) .or. &
           (amoeba .and. interaction == 'M' .and. field == 'P') .or. &
           (amoeba .and. interaction == 'S')) then
            ! Screening rules based on connectivity matrix: 
            ! 1. all the ones of non-amoeba FF, 
            ! 2. the static to IPD polarization field
            ! 3. the static-static interaction of AMOEBA
            do ineigh=1,4
                ! Look if j is at distance ineigh from i
                if(any(top%conn(ineigh)%ci(top%conn(ineigh)%ri(i_mm): &
                                           top%conn(ineigh)%ri(i_mm+1)-1) == j_mm)) then
                   
                    if(interaction == 'M' .and. amoeba) then
                        ! Amoeba uses two different screening rules for atoms 
                        ! of the same group and for atoms of different group
                        if(eel%mmat_polgrp(i_mm) == eel%mmat_polgrp(j_mm)) then
                            scalf = myscale_intra(ineigh)
                        else
                            scalf = myscale(ineigh)
                        end if
                    else
                        ! In any other case just use the scale factor for this
                        ! distance
                        scalf = myscale(ineigh)
                    end if
                    
                    ! Exit the loop
                    exit 

                end if
            end do
        else if((amoeba .and. interaction == 'M' .and. field == 'D') .or. &
                (amoeba .and. interaction == 'P')) then
            ! Screening rules based on polarization groups: 
            ! 1. the static to IPD direct field
            ! 2. the IPD to IPD interaction of AMOEBA
            pg_i = eel%mmat_polgrp(i_mm)

            outer: do ineigh=1,4
                do igrp=eel%pg_conn(ineigh)%ri(pg_i), &
                        eel%pg_conn(ineigh)%ri(pg_i+1)-1
                    ! Indexes of groups at distance ineigh from group pg_i
                    grp = eel%pg_conn(ineigh)%ci(igrp)
                    
                    if(any(eel%polgrp_mmat%ci(eel%polgrp_mmat%ri(grp): &
                                              eel%polgrp_mmat%ri(grp+1)-1) == j_mm)) then
                        ! If atom j is in a group at distance ineigh from the 
                        ! one of atom i, the field is scaled according to dscale
                        scalf = myscale(ineigh)
                        exit outer
                    end if
                end do
            end do outer
        else 
            ! Unexpected error
            call fatal_error("Unexpected combination of parameter for screening_rules")
        end if

    end function screening_rules