pitors_potential Subroutine

public subroutine pitors_potential(bds, V)

Uses

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

Compute pi-torsion terms of the potential.
This potential is defined on a -system, and uses the coordinates of six atoms A...F the central "double" bond is A-B, then C and D are connected to A while E and F are connected to B. So two plane ACD and BEF are defined. The potential is computed using the dihedral angle of the normal vector of those two planes, connected by segment A-B ().
The formula used is:

Arguments

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

pi-torsion potential, result will be added to V

Called by

proc~~pitors_potential~~CalledByGraph proc~pitors_potential pitors_potential proc~ommp_get_pitors_energy ommp_get_pitors_energy proc~ommp_get_pitors_energy->proc~pitors_potential proc~ommp_get_full_bnd_energy ommp_get_full_bnd_energy proc~ommp_get_full_bnd_energy->proc~pitors_potential proc~c_ommp_get_pitors_energy C_ommp_get_pitors_energy proc~c_ommp_get_pitors_energy->proc~ommp_get_pitors_energy program~test_si_potential test_SI_potential program~test_si_potential->proc~ommp_get_pitors_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_bnd_energy C_ommp_get_full_bnd_energy proc~c_ommp_get_full_bnd_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

pitors_potential

Source Code

    subroutine pitors_potential(bds, V)
        !! Compute pi-torsion terms of the potential.  
        !! This potential is defined on a \(\pi\)-system, and uses the 
        !! coordinates of six atoms A...F the central "double" bond is A-B, then
        !! C and D are connected to A while E and F are connected to B. So two
        !! plane ACD and BEF are defined. The potential is computed using the 
        !! dihedral angle of the normal vector of those two planes, connected 
        !! by segment A-B (\(\theta\)).  
        !! The formula used is:
        !! \[U_{\pi-torsion} = \sum_i k_i \large(1 + cos(2\theta-\pi) \large)\]

        use mod_constants, only : pi

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! pi-torsion potential, result will be added to V
        real(rp), dimension(3) :: a, b, c, d, e, f, u, t, cd, plv1, plv2, pln1, pln2
        real(rp) :: thet, costhet
        integer(ip) :: i

        if(.not. bds%use_pitors) return
        
        !$omp parallel do default(shared) reduction(+:V) &
        !$omp private(i,a,b,c,d,e,f,plv1,plv2,pln1,pln2,t,u,cd,thet,costhet)
        do i=1, bds%npitors
            !
            !  2(c)        5(e)         a => 1
            !   \         /             b => 4
            !    1(a) -- 4(b)  
            !   /         \
            !  3(d)        6(f)
            
            ! Atoms that defines the two planes
            a = bds%top%cmm(:,bds%pitorsat(1,i))
            c = bds%top%cmm(:,bds%pitorsat(2,i))
            d = bds%top%cmm(:,bds%pitorsat(3,i))

            b = bds%top%cmm(:,bds%pitorsat(4,i))
            e = bds%top%cmm(:,bds%pitorsat(5,i))
            f = bds%top%cmm(:,bds%pitorsat(6,i))


            ! Compute the normal vector of the first plane
            plv1 = d - b
            plv2 = c - b
            pln1(1) = plv1(2)*plv2(3) - plv1(3)*plv2(2)
            pln1(2) = plv1(3)*plv2(1) - plv1(1)*plv2(3)
            pln1(3) = plv1(1)*plv2(2) - plv1(2)*plv2(1)

            ! Compute the normal vector of the second plane
            plv1 = f - a
            plv2 = e - a
            pln2(1) = plv1(2)*plv2(3) - plv1(3)*plv2(2)
            pln2(2) = plv1(3)*plv2(1) - plv1(1)*plv2(3)
            pln2(3) = plv1(1)*plv2(2) - plv1(2)*plv2(1)

            cd = b - a

            t(1) = pln1(2)*cd(3) - pln1(3)*cd(2)
            t(2) = pln1(3)*cd(1) - pln1(1)*cd(3)
            t(3) = pln1(1)*cd(2) - pln1(2)*cd(1)
            t = t / norm2(t)
            
            u(1) = cd(2)*pln2(3) - cd(3)*pln2(2)
            u(2) = cd(3)*pln2(1) - cd(1)*pln2(3)
            u(3) = cd(1)*pln2(2) - cd(2)*pln2(1)
            u = u / norm2(u)
            
            costhet = dot_product(u,t)
            
            thet = acos(costhet)
            
            V = V +  bds%kpitors(i) * (1 + cos(2.0*thet-pi))
        end do

    end subroutine pitors_potential