mod_jacobian_mat.F90 Source File

module mod_jacobian_mat
    use mod_memory, only: rp, ip
    implicit none
    private

    public :: Rij_jacobian
    public :: simple_angle_jacobian, inplane_angle_jacobian, &
              torsion_angle_jacobian, opb_angle_jacobian, &
              pitors_angle_jacobian

    contains
        
    pure subroutine Rij_jacobian(ci, cj, Rij, J_i, J_j)
        !! Compute the Jacobian matrix of distance 
        !! Rij = sqrt((ci(_x_)-cj(_x_))**2 + (ci(_y_)-cj(_y_))**2 + (ci(_z_)-cj(_z_))**2)
        !! Derivatives wrt ci(:) are saved in J_i and wrt cj(:) in J_j; the
        !! distance between the two points is also provided in output in Rij.

        implicit none

        real(rp), intent(in) :: ci(3), cj(3)
        real(rp), intent(out) :: Rij, J_i(3), J_j(3)
        integer(ip) :: i

        Rij = norm2(ci-cj)
        do i=1, 3
            J_i(i) = (ci(i) - cj(i))/Rij
            J_j(i) = -J_i(i)
        end do
    end subroutine

    pure subroutine simple_angle_jacobian(ca, cb, cc, thet, J_a, J_b, J_c)
        implicit none

        real(rp), intent(in), dimension(3) :: ca, cb, cc
        !! Coordinates of the atoms defining the angle
        real(rp), intent(out), dimension(3) :: J_a, J_b, J_c
        !! The Jacobian components on atoms a, b and c respectively
        real(rp), intent(out) :: thet
        !! The angle (in rad) defined by ca-cb-cc

        real(rp), dimension(3) :: dr1, dr2
        real(rp) :: l1, l2, dr1_d_dr2, acosd, cost
        integer(ip) :: i

        dr1 = ca - cb
        dr2 = cc - cb
        dr1_d_dr2 = dot_product(dr1, dr2)
        l1 = norm2(dr1)
        l2 = norm2(dr2)
        cost = dr1_d_dr2/(l1*l2)
        thet = acos(cost)
        acosd = 1.0 / sqrt(1.0 - cost**2)

        do i=1,3
            J_a(i) = -acosd * (dr2(i) * l1 - dr1(i) * dr1_d_dr2 / l1) / (l1**2*l2)
            J_c(i) = -acosd * (dr1(i) * l2 - dr2(i) * dr1_d_dr2 / l2) / (l1*l2**2)
            J_b(i) = -(J_a(i) + J_c(i))
        end do
    end subroutine
    
    pure subroutine inplane_angle_jacobian(ca, cb, cc, cx, thet, &
                                           J_a, J_b, J_c, J_x)
        !! Computes the Jacobian matrix for the inplane angle definition.
        !! It computes the Jacobian for the normal angle using the projected 
        !! point (R) as central point. Then projects 
        !! \(J_r = \frac{\partial \theta}{\partial \vec{R}} \) onto A, B, C, and X 
        !! (auxiliary point). The projection is done computing the 3x3 matrices
        !! of partial derivative of \(\vec{R}\) wrt any actual point and using
        !! them to project \(J_r\).
        !! \[\frac{\partial \vec{R}}{\partial \vec{A}} = 
        !!      \begin{bmatrix} 
        !!          \frac{\partial \vec{R}_x}{\partial \vec{A}_x} & 
        !!          \frac{\partial \vec{R}_y}{\partial \vec{A}_x} & 
        !!          \frac{\partial \vec{R}_z}{\partial \vec{A}_x} \\ 
        !!          \frac{\partial \vec{R}_x}{\partial \vec{A}_y} & 
        !!          \frac{\partial \vec{R}_y}{\partial \vec{A}_y} & 
        !!          \frac{\partial \vec{R}_z}{\partial \vec{A}_y} \\ 
        !!          \frac{\partial \vec{R}_x}{\partial \vec{A}_z} & 
        !!          \frac{\partial \vec{R}_y}{\partial \vec{A}_z} & 
        !!          \frac{\partial \vec{R}_z}{\partial \vec{A}_z} \\ 
        !!      \end{bmatrix} \]
        !! Those matrices are computed using the chain rule.
        !!
        !! Exemple:
        !! \[\vec{V} = \vec{B} - \vec{X}\]
        !! \[\vec{P} = (\vec{A} - \vec{X}) \times (\vec{C} - \vec{X})\]
        !! \[\hat{P} = \frac{\vec{P}}{||\vec{P}||}\]
        !! \[\vec{R} = \vec{B} - (\vec{V}\cdot\hat{P})\hat{P}\]
        !! \[J_a = \frac{\partial \theta}{\partial \vec{A}} 
        !! = \left( \frac{\partial \theta}{\partial \vec{A}} \right)_\vec{R} + 
        !!   \frac{\partial \vec{R}}{\partial \vec{A}} \times
        !!   J_r
        !! = \left( \frac{\partial \theta}{\partial \vec{A}} \right)_\vec{R} +
        !!   - \frac{\partial \vec{P}}{\partial \vec{A}} \times
        !!     \frac{\partial \hat{P}}{\partial \vec{P}} \times 
        !!     \frac{\partial (\vec{V} \cdot \hat{P})}{\partial \hat{P}} \times
        !!     J_r
        !! \]
        use mod_utils, only: cross_product, vec_skw, versor_der

        implicit none

        real(rp), intent(in), dimension(3) :: ca, cb, cc, cx
        !! Coordinates of the atoms defining the angle
        real(rp), intent(out), dimension(3) :: J_a, J_b, J_c, J_x
        !! The Jacobian components on atoms a, b and c respectively
        real(rp), intent(out) :: thet
        !! The angle (in rad) defined by ca-cb-cc

        real(rp), dimension(3) :: cr, cv, cp, cpp, cq, cs, J_r
        real(rp), dimension(3,3) :: drda, drdc, drdx, dkpdp, dppda, &
                                    dpdpp, dpda, dpdc, dppdc, dpdx, dppdx, &
                                    dkpdv
        real(rp) :: k
        integer(ip) :: i

        cq = ca - cx
        cs = cc - cx
        cv = cb - cx

        cpp = cross_product(cq, cs)
        cp = cpp / norm2(cpp)

        k = dot_product(cv, cp)
        cr = cb - k * cp
        
        call simple_angle_jacobian(ca, cr, cc, thet, J_a, J_r, J_c)

        do i=1,3
            dkpdp(i,:) = cv(i) * cp
            dkpdp(i,i) = dkpdp(i,i) + k
        end do

        dpdpp = versor_der(cpp)

        do i=1,3
            dkpdv(i,:) = cp(i) * cp
        end do

        dppda = -vec_skw(cs)
        dppdc = vec_skw(cq)
       
        dppdx = -vec_skw(cs-cq)

        dpda = matmul(dppda, dpdpp)
        drda = matmul(dpda, dkpdp)

        dpdc = matmul(dppdc, dpdpp)
        drdc = matmul(dpdc, dkpdp)

        dpdx = matmul(dppdx, dpdpp)
        drdx = -(matmul(dpdx, dkpdp) - dkpdv)

        J_a = J_a + matmul(drda, J_r)
        J_x = matmul(drdx, J_r)
        J_c = J_c + matmul(drdc, J_r)
        J_b = -(J_a+J_c+J_x)
        
    end subroutine

    subroutine torsion_angle_jacobian(ca, cb, cc, cd, thet, &
                                           J_a, J_b, J_c, J_d)
        !! Computes the Jacobian matrix for torsion angle defined by points 
        !! \(\vec{A}\), \(\vec{B}\), \(\vec{C}\) and \(\vec{D}\) (connected in
        !! this order). The angle is defined as follow:
        !! \[ \vec{U} = (\vec{B} - \vec{C}) \times (\vec{D} - \vec{C}) \\
        !!    \vec{T} = (\vec{B} - \vec{A}) \times (\vec{B} - \vec{C}) \\
        !!    cos(\theta) = \vec{U} \cdot \vec{T}
        !! \]
        !! Applying the chain rule:
        !! \[J_a = \frac{\partial \theta}{\partial \vec{A}} 
        !!       = -\frac{1}{\sqrt{1-\theta^2}} 
        !!          \frac{\partial cos(\theta)}{\partial \vec{A}} 
        !!       = -\frac{1}{\sqrt{1-\theta^2}}
        !!          \frac{\partial \vec{U} \cdot \vec{T}}{\partial \vec{A}} \\
        !!       = -\frac{1}{\sqrt{1-\theta^2}}
        !!          \frac{\partial (\vec{B} - \vec{A})}{\partial \vec{A}} \times
        !!          \frac{\partial \vec{T}}{\partial (\vec{B} - \vec{A})} \times
        !!          \frac{\partial \hat{T}}{\partial \vec{T}} \times \vec{U} \\
        !!       = -\frac{1}{\sqrt{1-\theta^2}}
        !!          (- \mathbb{I}_3) \times
        !!          skw(B-C) \times
        !!          \frac{\partial \hat{T}}{\partial \vec{T}} \times \vec{U}
        !! \]
        !! \[J_d = \frac{\partial \theta}{\partial \vec{D}} =
        !!        = -\frac{1}{\sqrt{1-\theta^2}}
        !!          \mathbb{I}_3 \times (-skw(B-C)) \times
        !!          \frac{\partial \hat{U}}{\partial \vec{U}} \times \vec{T}
        !! \]
        use mod_utils, only: cross_product, vec_skw, versor_der
        use mod_constants, only: pi, eps_rp

        implicit none

        real(rp), intent(inout), dimension(3) :: ca, cb, cc, cd
        !! Coordinates of the atoms defining the angle
        real(rp), intent(out), dimension(3) :: J_a, J_b, J_c, J_d
        !! The Jacobian components on atoms a, b and c and d respectively
        real(rp), intent(out) :: thet
        !! The torsion angle

        real(rp), dimension(3) :: a_b, c_d, c_b, t, u, ht, hu
        real(rp), dimension(3,3) :: dhudu, dhtdt, dtda, dudd, dhudd, dhtda, &
                                    dhudb, dhtdb, dudb, dtdb, &
                                    dhudc, dtdc, dudc
        real(rp) :: costhet, dacost, s

        a_b = cb - ca
        c_d = cd - cc
        c_b = cb - cc

        t = cross_product(a_b,c_b)
        ht = t / norm2(t)
        
        u = cross_product(c_b,c_d)
        hu = u / norm2(u)

        costhet = dot_product(hu,ht)
        if(costhet + 1.0 <= eps_rp) then
            thet = pi
            dacost = -1.0 / sin(costhet)
        else if(abs(costhet - 1.0) <= eps_rp) then
            thet = 0.0
            dacost = -1.0 / sin(costhet)
        else
            thet = acos(costhet)
            dacost = - 1.0/sqrt(1.0-costhet**2)
        end if

        s = sign(1.0_rp, -dot_product(a_b, u))
        thet = thet * s
        
        dacost = dacost * s
        
        dhtdt = versor_der(t)
        dhudu = versor_der(u)

        dtda = vec_skw(c_b)
        dudd = dtda 
        
        dhtda = matmul(dtda, dhtdt)
        J_a = -dacost * matmul(dhtda,hu) 


        dudb = vec_skw(c_d)
        dtdb = vec_skw(ca-cc)

        dhudb = matmul(dudb, dhudu)
        dhtdb = matmul(dtdb, dhtdt)
        J_b = dacost * (matmul(dhudb,ht) + matmul(dhtdb,hu))

        dtdc = vec_skw(a_b)
        dudc = vec_skw(cb - cd)

        dhtdb = matmul(dtdc, dhtdt)
        dhudc = matmul(dudc, dhudu) 
        J_c = dacost * (matmul(dhudc,ht) + matmul(dhtdb,hu))

        dhudd = matmul(dudd, dhudu)
        J_d = -dacost * matmul(dhudd,ht)

    end subroutine
    
    subroutine opb_angle_jacobian(ca, cb, cc, cd, thet, &
                                  J_a, J_b, J_c, J_d)
        use mod_utils, only: cross_product, vec_skw, versor_der
        use mod_constants, only: pi, eps_rp

        implicit none

        real(rp), intent(inout), dimension(3) :: ca, cb, cc, cd
        !! Coordinates of the atoms defining the angle
        real(rp), intent(out), dimension(3) :: J_a, J_b, J_c, J_d
        !! The Jacobian components on atoms a, b and c and d respectively
        real(rp), intent(out) :: thet
        !! The out-of-plane angle

        real(rp), dimension(3) :: a_b, a_c, v, p, hp, hv
        real(rp), dimension(3,3) :: dhpdp, dhvdv, &
                                    dhpdb, dpdb, dhpdc, dpdc, dvdd
        real(rp) :: costhet, dacost, thet0

        a_b = ca - cb
        a_c = ca - cc
        v = ca - cd

        p = cross_product(a_b,a_c)
        
        hp = p / norm2(p)
        hv = v / norm2(v)
        
        costhet = dot_product(hv, hp)
        if(costhet + 1.0 <= eps_rp) then
            thet0 = pi
        else 
            thet0 = acos(costhet)
        end if

        thet = abs(pi/2.0 - thet0)
        dacost = 1.0/sqrt(1.0-costhet**2)
        if(pi/2.0 - thet0 < 0.0) dacost = -dacost
        
        dhpdp = versor_der(p)
        dhvdv = versor_der(v)
        
        dpdb = vec_skw(a_c)
        dhpdb = matmul(dpdb,dhpdp)
        J_b = -dacost * matmul(dhpdb,hv)

        dpdc = vec_skw(a_b)
        dpdc = -dpdc
        dhpdc = matmul(dpdc, dhpdp)
        J_c = -dacost * matmul(dhpdc,hv) 

        dvdd = -dhvdv
        J_d = dacost * matmul(dvdd,hp)
        
        J_a = -(J_b + J_c +J_d)
    end subroutine

    subroutine pitors_angle_jacobian(ca, cb, cc, cd, ce, cf, thet, &
                                     J_a, J_b, J_c, J_d, J_e, J_f)
        use mod_utils, only: cross_product, vec_skw, versor_der

        implicit none

        real(rp), dimension(3), intent(inout) :: ca, cb, cc, cd, ce, cf
        real(rp), intent(out) :: thet
        real(rp), dimension(3), intent(out) :: J_a, J_b, J_c, J_d, J_e, J_f

        real(rp), dimension(3) :: d_b, c_b, f_a, e_a, s, p, r, t, ht, u, hu
        real(rp) :: costhet, dcosthet

        real(rp), dimension(3,3) :: skw_s, dhtdt, dhudu, &
                                    skw_c_b, skw_d_b, skw_f_a, skw_e_a, &
                                    skw_r, skw_p, dpdb, dtdb
        real(rp), dimension(3) :: dcostdp, dcostdr, dcostdt, dcostdu

        d_b = cd - cb
        c_b = cc - cb
        f_a = cf - ca
        e_a = ce - ca
        s = cb - ca
        p = cross_product(d_b, c_b)
        r = cross_product(f_a, e_a)

        t = cross_product(p, s)
        ht = t / norm2(t)
        u = cross_product(s, r)
        hu = u / norm2(u)

        costhet = dot_product(hu, ht)
        thet = acos(costhet)
        dcosthet = - 1.0 / sqrt(1.0 - costhet**2)
        
        skw_s = vec_skw(s)
        skw_r = vec_skw(r)
        skw_p = vec_skw(p)
        
        dhudu = versor_der(u)
        dhtdt = versor_der(t)

        dcostdt = matmul(dhtdt, hu)
        dcostdu = matmul(dhudu, ht)
        
        dcostdp = matmul(skw_s, dcostdt)
        dcostdr = -matmul(skw_s, dcostdu)

        skw_c_b = vec_skw(c_b)
        skw_d_b = vec_skw(d_b)
        skw_e_a = vec_skw(e_a)
        skw_f_a = vec_skw(f_a)

        dpdb = skw_d_b - skw_c_b
        dtdb = matmul(dpdb, skw_s) - skw_p

        J_b = dcosthet * (matmul(skw_r, dcostdu) + matmul(dtdb, dcostdt))
        J_c = dcosthet * (-matmul(skw_d_b, dcostdp))
        J_d = dcosthet * (matmul(skw_c_b, dcostdp))
        J_e = dcosthet * (-matmul(skw_f_a, dcostdr))
        J_f = dcosthet * (matmul(skw_e_a, dcostdr))

        J_a = -(J_b + J_c + J_d + J_e + J_f)
    end subroutine
end module mod_jacobian_mat