bond_potential Subroutine

public subroutine bond_potential(bds, V)

Uses

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

Compute the bond-stretching terms of the potential energy.
They are computed according to the general formula adopted in AMOEBA Force Field:

Arguments

Type IntentOptional Attributes Name
type(ommp_bonded_type), intent(in) :: bds
real(kind=rp), intent(inout) :: V

Bond potential, result will be added to V

Called by

proc~~bond_potential~~CalledByGraph proc~bond_potential bond_potential proc~ommp_get_bond_energy ommp_get_bond_energy proc~ommp_get_bond_energy->proc~bond_potential proc~ommp_get_full_bnd_energy ommp_get_full_bnd_energy proc~ommp_get_full_bnd_energy->proc~bond_potential proc~c_ommp_get_bond_energy C_ommp_get_bond_energy proc~c_ommp_get_bond_energy->proc~ommp_get_bond_energy proc~c_ommp_get_full_bnd_energy C_ommp_get_full_bnd_energy proc~c_ommp_get_full_bnd_energy->proc~ommp_get_full_bnd_energy proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_bnd_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

bond_potential

Source Code

    subroutine bond_potential(bds, V)
        !! Compute the bond-stretching terms of the potential energy.  
        !! They are computed according to the general formula adopted in AMOEBA
        !! Force Field:
        !! \[U_{bond} = \sum_i k_i \Delta l_i^2 \large(1 + k^{(3)}\Delta l_i + 
        !! k^{(4)}\Delta l_i^2 \large)\]
        !! \[\Delta l_i = l_i - l^{(eq)}_i\]

        use mod_constants, only : eps_rp

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! Bond potential, result will be added to V

        integer :: i
        logical(lp) :: use_cubic, use_quartic
        real(rp) :: dr(3), l, dl, dl2

        use_cubic = (abs(bds%bond_cubic) > eps_rp)
        use_quartic = (abs(bds%bond_quartic) > eps_rp)
        
        if(.not. bds%use_bond) return

        if(.not. use_cubic .and. .not. use_quartic) then
            ! This is just a regular harmonic potential
            !$omp parallel do default(shared) schedule(static) & 
            !$omp private(i,dr,l,dl) reduction(+:v) 
            do i=1, bds%nbond
                dr = bds%top%cmm(:,bds%bondat(1,i)) - &
                     bds%top%cmm(:,bds%bondat(2,i))
                l = sqrt(dot_product(dr, dr))
                dl = l - bds%l0bond(i)
                
                V = V + bds%kbond(i) * dl * dl
            end do
        else
            !$omp parallel do default(shared) schedule(static) & 
            !$omp private(i,dr,l,dl,dl2) reduction(+:v) 
            do i=1, bds%nbond
                dr = bds%top%cmm(:,bds%bondat(1,i)) - &
                     bds%top%cmm(:,bds%bondat(2,i))
                l = sqrt(dot_product(dr, dr))
                dl = l - bds%l0bond(i)
                dl2 = dl * dl

                V = V + bds%kbond(i)*dl2 * &
                    (1.0_rp + bds%bond_cubic*dl + bds%bond_quartic*dl2)
            end do
        end if
        
    end subroutine bond_potential