pitors_potential Subroutine

    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