strbnd_potential Subroutine

public subroutine strbnd_potential(bds, V)

Compute the stretch-bend cross term potential.
Those terms are computed according the following formula: where is the angle delimited by the bond and .
The force constants and are explicitely defined in the FF, while the equilibrium values are the same as for stretching and bending terms.

Arguments

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

Stretch-bend cross term potential, result will be added to V

Called by

proc~~strbnd_potential~~CalledByGraph proc~strbnd_potential strbnd_potential proc~ommp_get_strbnd_energy ommp_get_strbnd_energy proc~ommp_get_strbnd_energy->proc~strbnd_potential proc~ommp_get_full_bnd_energy ommp_get_full_bnd_energy proc~ommp_get_full_bnd_energy->proc~strbnd_potential proc~c_ommp_get_strbnd_energy C_ommp_get_strbnd_energy proc~c_ommp_get_strbnd_energy->proc~ommp_get_strbnd_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

strbnd_potential

Source Code

    subroutine strbnd_potential(bds, V)
        !! Compute the stretch-bend cross term potential.   
        !! Those terms are computed according the following formula:
        !! \[U_{bond/angle} = (k_i \Delta l_i + k_j \Delta l_j) 
        !! \Delta \theta_{ij} \]
        !! where \(\theta_{ij}\) is the angle delimited by the bond \(i\) and 
        !! \(j\).   
        !! The force constants \(k_i\) and \(k_j\) are explicitely defined in
        !! the FF, while the equilibrium values are the same as for stretching
        !! and bending terms.

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! Stretch-bend cross term potential, result will be added to V

        integer(ip) :: i
        real(rp) :: d_l1, d_l2, d_thet, dr1(3), dr2(3), l1, l2, thet
        
        if(.not. bds%use_strbnd) return

        !$omp parallel do default(shared) reduction(+:V) &
        !$omp private(i,dr1,l1,l2,d_l1,d_l2,dr2,thet,d_thet)
        do i=1, bds%nstrbnd
            dr1 = bds%top%cmm(:, bds%strbndat(2,i)) - &
                  bds%top%cmm(:, bds%strbndat(1,i))
            l1 = norm2(dr1)
            d_l1 = l1 - bds%strbndl10(i)
            
            dr2 = bds%top%cmm(:, bds%strbndat(2,i)) - &
                  bds%top%cmm(:, bds%strbndat(3,i))
            l2 = norm2(dr2)
            d_l2 = l2 - bds%strbndl20(i)

            thet = acos(dot_product(dr1, dr2)/(l1*l2))
            d_thet = thet - bds%strbndthet0(i) 
            
            V = V + (d_l1*bds%strbndk1(i) + d_l2*bds%strbndk2(i)) * d_thet
        end do
    end subroutine strbnd_potential