opb_potential Subroutine

public subroutine opb_potential(bds, V)

Uses

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

Computes the out-of-plane bending potential.
With Allinger formula: similarly to in plane angles, here we are considering a trigonal center, where D is the central atom and A, B, C are connected to D. Allinger formula consider the angle between vector and the normal vector of plane ABC, using as implicit equilibrium value. The formula for this potential term is:

Arguments

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

out-of-plane potential, result will be added to V

Called by

proc~~opb_potential~~CalledByGraph proc~opb_potential opb_potential proc~ommp_get_opb_energy ommp_get_opb_energy proc~ommp_get_opb_energy->proc~opb_potential proc~ommp_get_full_bnd_energy ommp_get_full_bnd_energy proc~ommp_get_full_bnd_energy->proc~opb_potential proc~c_ommp_get_opb_energy C_ommp_get_opb_energy proc~c_ommp_get_opb_energy->proc~ommp_get_opb_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

opb_potential

Source Code

    subroutine opb_potential(bds, V)
        !! Computes the out-of-plane bending potential.  
        !! With Allinger formula: similarly to in plane angles, here we are 
        !! considering a trigonal center, where D is the central atom and 
        !! A, B, C are connected to D. Allinger formula consider the angle 
        !! between vector \(\vec{AD}\) and the normal vector of plane ABC, 
        !! using \(\frac{\pi}{2}\) as implicit equilibrium value. The formula
        !! for this potential term is:
        !! \[U_{out-of-plane} = \sum_i k_i \chi_i^2 \large(1 + 
        !! \sum_{j=1}^4 k^{(j+2)} \chi_i^j \large) \]

        use mod_constants, only : pi

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! out-of-plane potential, result will be added to V
        real(rp), dimension(3) :: a, b, c, d, plv1, plv2, pln, vad
        real(rp) :: lpln, lvad, thet, thet2, thet3, thet4
        integer(ip) :: i

        if(.not. bds%use_opb) return
        
        !$omp parallel do default(shared) reduction(+:V) &
        !$omp private(i,a,b,c,d,plv1,plv2,pln,lpln,vad,lvad,thet,thet2,thet3,thet4)
        do i=1, bds%nopb
            ! A* -- D -- C
            !       |
            !       B 
            a = bds%top%cmm(:,bds%opbat(2,i))
            d = bds%top%cmm(:,bds%opbat(1,i))
            c = bds%top%cmm(:,bds%opbat(3,i))
            b = bds%top%cmm(:,bds%opbat(4,i))

            ! Compute the normal vector of the plane
            plv1 = a - b
            plv2 = a - c
            pln(1) = plv1(2)*plv2(3) - plv1(3)*plv2(2)
            pln(2) = plv1(3)*plv2(1) - plv1(1)*plv2(3)
            pln(3) = plv1(1)*plv2(2) - plv1(2)*plv2(1)
            lpln = norm2(pln)

            ! Vector from A to D
            vad = a - d
            lvad = norm2(vad)

            thet = abs(pi/2.0 - acos(dot_product(vad, pln)/(lvad*lpln)))
            thet2 = thet*thet
            thet3 = thet2*thet
            thet4 = thet3*thet
            V = V +  bds%kopb(i) * thet2 * (1 + bds%opb_cubic*thet &
                + bds%opb_quartic*thet2 + bds%opb_pentic*thet3 &
                + bds%opb_sextic*thet4)
        end do
    end subroutine opb_potential