Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(ommp_bonded_type), | intent(in) | :: | bds | |||
real(kind=rp), | intent(inout) | :: | grad(3,bds%top%mm_atoms) |
Gradients of bond stretching terms of potential energy |
subroutine strbnd_geomgrad(bds, grad)
use mod_jacobian_mat, only: Rij_jacobian, simple_angle_jacobian
implicit none
type(ommp_bonded_type), intent(in) :: bds
! Bonded potential data structure
real(rp), intent(inout) :: grad(3,bds%top%mm_atoms)
!! Gradients of bond stretching terms of potential energy
integer(ip) :: i, ia, ib, ic
real(rp) :: d_l1, d_l2, d_thet, l1, l2, thet, g1, g2, g3
real(rp), dimension(3) :: a, b, c, &
J1_a, J1_b, &
J2_b, J2_c, &
J3_a, J3_b, J3_c
logical :: sk_a, sk_b, sk_c
if(.not. bds%use_strbnd) return
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,ia,ib,ic,sk_a,sk_b,sk_c,a,b,c,l1,l2,d_l1,d_l2,thet,d_thet) &
!$omp private(J1_a,J1_b,J2_b,J2_c,J3_a,J3_b,J3_c,g1,g2,g3)
do i=1, bds%nstrbnd
ia = bds%strbndat(1,i)
ib = bds%strbndat(2,i)
ic = bds%strbndat(3,i)
if(bds%top%use_frozen) then
sk_a = bds%top%frozen(ia)
sk_b = bds%top%frozen(ib)
sk_c = bds%top%frozen(ic)
if(sk_a .and. sk_b .and. sk_c) cycle
else
sk_a = .false.
sk_b = .false.
sk_c = .false.
end if
a = bds%top%cmm(:, ia)
b = bds%top%cmm(:, ib)
c = bds%top%cmm(:, ic)
call Rij_jacobian(a, b, l1, J1_a, J1_b)
call Rij_jacobian(b, c, l2, J2_b, J2_c)
call simple_angle_jacobian(a, b, c, thet, J3_a, J3_b, J3_c)
d_l1 = l1 - bds%strbndl10(i)
d_l2 = l2 - bds%strbndl20(i)
d_thet = thet - bds%strbndthet0(i)
g1 = bds%strbndk1(i) * d_thet
g2 = bds%strbndk2(i) * d_thet
g3 = bds%strbndk1(i) * d_l1 + bds%strbndk2(i) * d_l2
if(.not. sk_a) then
!$omp atomic update
grad(1,ia) = grad(1,ia) + J1_a(1) * g1 + J3_a(1) * g3
!$omp atomic update
grad(2,ia) = grad(2,ia) + J1_a(2) * g1 + J3_a(2) * g3
!$omp atomic update
grad(3,ia) = grad(3,ia) + J1_a(3) * g1 + J3_a(3) * g3
end if
if(.not. sk_b) then
!$omp atomic update
grad(1,ib) = grad(1,ib) + J1_b(1) * g1 + J2_b(1) * g2 + J3_b(1) * g3
!$omp atomic update
grad(2,ib) = grad(2,ib) + J1_b(2) * g1 + J2_b(2) * g2 + J3_b(2) * g3
!$omp atomic update
grad(3,ib) = grad(3,ib) + J1_b(3) * g1 + J2_b(3) * g2 + J3_b(3) * g3
end if
if(.not. sk_c) then
!$omp atomic update
grad(1,ic) = grad(1,ic) + J2_c(1) * g2 + J3_c(1) * g3
!$omp atomic update
grad(2,ic) = grad(2,ic) + J2_c(2) * g2 + J3_c(2) * g3
!$omp atomic update
grad(3,ic) = grad(3,ic) + J2_c(3) * g2 + J3_c(3) * g3
end if
end do
end subroutine strbnd_geomgrad