| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_bonded_type), | intent(in) | :: | bds | |||
| real(kind=rp), | intent(inout) | :: | grad(3,bds%top%mm_atoms) |
improper torsion potential, result will be added to V |
subroutine strtor_geomgrad(bds, grad)
use mod_jacobian_mat, only: Rij_jacobian, torsion_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)
!! improper torsion potential, result will be added to V
real(rp) :: thet, gt(3), dihef(3), dr1, dr2, dr3, r1, r2, r3, &
Jt_a(3), Jt_b(3), Jt_c(3), Jt_d(3), &
Jb1_a(3), Jb1_b(3), &
Jb2_a(3), Jb2_b(3), &
Jb3_a(3), Jb3_b(3)
integer(ip) :: i, j, k, ib1, ib2, ib3, &
it_a, it_b, it_c, it_d, &
ib1_a, ib1_b, &
ib2_a, ib2_b, &
ib3_a, ib3_b
logical :: sk_ta, sk_tb, sk_tc, sk_td, &
sk_1a, sk_1b, &
sk_2a, sk_2b, &
sk_3a, sk_3b
if (.not. bds%use_strtor) return
!$omp parallel do default(shared) &
!$omp private(thet, gt, dihef, dr1, dr2, dr3, r1, r2, r3) &
!$omp private(Jt_a, Jt_b, Jt_c, Jt_d, Jb1_a, Jb1_b, Jb2_a, Jb2_b) &
!$omp private(Jb3_a, Jb3_b, i, j, k, ib1, ib2, ib3, it_a, it_b, it_c, it_d) &
!$omp private(ib1_a, ib1_b, ib2_a, ib2_b, ib3_a, ib3_b, sk_ta, sk_tb, sk_tc, sk_td) &
!$omp private(sk_1a, sk_1b, sk_2a, sk_2b, sk_3a, sk_3b)
do i = 1, bds%nstrtor
! Atoms that define the dihedral angle
it_a = bds%strtorat(1, i)
it_b = bds%strtorat(2, i)
it_c = bds%strtorat(3, i)
it_d = bds%strtorat(4, i)
ib1 = bds%strtor_b(1, i)
ib1_a = bds%bondat(1, ib1)
ib1_b = bds%bondat(2, ib1)
ib2 = bds%strtor_b(2, i)
ib2_a = bds%bondat(1, ib2)
ib2_b = bds%bondat(2, ib2)
ib3 = bds%strtor_b(3, i)
ib3_a = bds%bondat(1, ib3)
ib3_b = bds%bondat(2, ib3)
if (bds%top%use_frozen) then
sk_ta = bds%top%frozen(it_a)
sk_tb = bds%top%frozen(it_b)
sk_tc = bds%top%frozen(it_c)
sk_td = bds%top%frozen(it_d)
sk_1a = bds%top%frozen(ib1_a)
sk_1b = bds%top%frozen(ib1_b)
sk_2a = bds%top%frozen(ib2_a)
sk_2b = bds%top%frozen(ib2_b)
sk_3a = bds%top%frozen(ib3_a)
sk_3b = bds%top%frozen(ib3_b)
if (sk_ta .and. sk_tb .and. sk_tc .and. sk_td .and. &
sk_1a .and. sk_1b .and. &
sk_2a .and. sk_2b .and. &
sk_3a .and. sk_3b) cycle
else
sk_ta = .false.
sk_tb = .false.
sk_tc = .false.
sk_td = .false.
sk_1a = .false.
sk_1b = .false.
sk_2a = .false.
sk_2b = .false.
sk_3a = .false.
sk_3b = .false.
end if
call torsion_angle_jacobian(bds%top%cmm(:, it_a), &
bds%top%cmm(:, it_b), &
bds%top%cmm(:, it_c), &
bds%top%cmm(:, it_d), &
thet, Jt_a, Jt_b, Jt_c, Jt_d)
do j = 1, 3
gt(j) = -real(j) * sin(j * thet + bds%torsphase(j, bds%angtor_t(i)))
dihef(j) = 1.0 + cos(j * thet + bds%torsphase(j, bds%angtor_t(i)))
end do
call Rij_jacobian(bds%top%cmm(:, ib1_a), &
bds%top%cmm(:, ib1_b), &
r1, Jb1_a, Jb1_b)
dr1 = r1 - bds%l0bond(ib1)
call Rij_jacobian(bds%top%cmm(:, ib2_a), &
bds%top%cmm(:, ib2_b), &
r2, Jb2_a, Jb2_b)
dr2 = r2 - bds%l0bond(ib2)
call Rij_jacobian(bds%top%cmm(:, ib3_a), &
bds%top%cmm(:, ib3_b), &
r3, Jb3_a, Jb3_b)
dr3 = r3 - bds%l0bond(ib3)
do k = 1, 3
if (.not. sk_ta) then
!$omp atomic update
grad(1, it_a) = grad(1, it_a) + bds%strtork(k, i) * dr1 * gt(k) * Jt_a(1)
!$omp atomic update
grad(2, it_a) = grad(2, it_a) + bds%strtork(k, i) * dr1 * gt(k) * Jt_a(2)
!$omp atomic update
grad(3, it_a) = grad(3, it_a) + bds%strtork(k, i) * dr1 * gt(k) * Jt_a(3)
end if
if (.not. sk_tb) then
!$omp atomic update
grad(1, it_b) = grad(1, it_b) + bds%strtork(k, i) * dr1 * gt(k) * Jt_b(1)
!$omp atomic update
grad(2, it_b) = grad(2, it_b) + bds%strtork(k, i) * dr1 * gt(k) * Jt_b(2)
!$omp atomic update
grad(3, it_b) = grad(3, it_b) + bds%strtork(k, i) * dr1 * gt(k) * Jt_b(3)
end if
if (.not. sk_tc) then
!$omp atomic update
grad(1, it_c) = grad(1, it_c) + bds%strtork(k, i) * dr1 * gt(k) * Jt_c(1)
!$omp atomic update
grad(2, it_c) = grad(2, it_c) + bds%strtork(k, i) * dr1 * gt(k) * Jt_c(2)
!$omp atomic update
grad(3, it_c) = grad(3, it_c) + bds%strtork(k, i) * dr1 * gt(k) * Jt_c(3)
end if
if (.not. sk_td) then
!$omp atomic update
grad(1, it_d) = grad(1, it_d) + bds%strtork(k, i) * dr1 * gt(k) * Jt_d(1)
!$omp atomic update
grad(2, it_d) = grad(2, it_d) + bds%strtork(k, i) * dr1 * gt(k) * Jt_d(2)
!$omp atomic update
grad(3, it_d) = grad(3, it_d) + bds%strtork(k, i) * dr1 * gt(k) * Jt_d(3)
end if
if (.not. sk_1a) then
!$omp atomic update
grad(1, ib1_a) = grad(1, ib1_a) + bds%strtork(k, i) * dihef(k) * Jb1_a(1)
!$omp atomic update
grad(2, ib1_a) = grad(2, ib1_a) + bds%strtork(k, i) * dihef(k) * Jb1_a(2)
!$omp atomic update
grad(3, ib1_a) = grad(3, ib1_a) + bds%strtork(k, i) * dihef(k) * Jb1_a(3)
end if
if (.not. sk_1b) then
!$omp atomic update
grad(1, ib1_b) = grad(1, ib1_b) + bds%strtork(k, i) * dihef(k) * Jb1_b(1)
!$omp atomic update
grad(2, ib1_b) = grad(2, ib1_b) + bds%strtork(k, i) * dihef(k) * Jb1_b(2)
!$omp atomic update
grad(3, ib1_b) = grad(3, ib1_b) + bds%strtork(k, i) * dihef(k) * Jb1_b(3)
end if
if (.not. sk_ta) then
!$omp atomic update
grad(1, it_a) = grad(1, it_a) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_a(1)
!$omp atomic update
grad(2, it_a) = grad(2, it_a) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_a(2)
!$omp atomic update
grad(3, it_a) = grad(3, it_a) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_a(3)
end if
if (.not. sk_tb) then
!$omp atomic update
grad(1, it_b) = grad(1, it_b) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_b(1)
!$omp atomic update
grad(2, it_b) = grad(2, it_b) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_b(2)
!$omp atomic update
grad(3, it_b) = grad(3, it_b) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_b(3)
end if
if (.not. sk_tc) then
!$omp atomic update
grad(1, it_c) = grad(1, it_c) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_c(1)
!$omp atomic update
grad(2, it_c) = grad(2, it_c) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_c(2)
!$omp atomic update
grad(3, it_c) = grad(3, it_c) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_c(3)
end if
if (.not. sk_td) then
!$omp atomic update
grad(1, it_d) = grad(1, it_d) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_d(1)
!$omp atomic update
grad(2, it_d) = grad(2, it_d) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_d(2)
!$omp atomic update
grad(3, it_d) = grad(3, it_d) + bds%strtork(3 + k, i) * dr2 * gt(k) * Jt_d(3)
end if
if (.not. sk_2a) then
!$omp atomic update
grad(1, ib2_a) = grad(1, ib2_a) + bds%strtork(3 + k, i) * dihef(k) * Jb2_a(1)
!$omp atomic update
grad(2, ib2_a) = grad(2, ib2_a) + bds%strtork(3 + k, i) * dihef(k) * Jb2_a(2)
!$omp atomic update
grad(3, ib2_a) = grad(3, ib2_a) + bds%strtork(3 + k, i) * dihef(k) * Jb2_a(3)
end if
if (.not. sk_2b) then
!$omp atomic update
grad(1, ib2_b) = grad(1, ib2_b) + bds%strtork(3 + k, i) * dihef(k) * Jb2_b(1)
!$omp atomic update
grad(2, ib2_b) = grad(2, ib2_b) + bds%strtork(3 + k, i) * dihef(k) * Jb2_b(2)
!$omp atomic update
grad(3, ib2_b) = grad(3, ib2_b) + bds%strtork(3 + k, i) * dihef(k) * Jb2_b(3)
end if
if (.not. sk_ta) then
!$omp atomic update
grad(1, it_a) = grad(1, it_a) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_a(1)
!$omp atomic update
grad(2, it_a) = grad(2, it_a) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_a(2)
!$omp atomic update
grad(3, it_a) = grad(3, it_a) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_a(3)
end if
if (.not. sk_tb) then
!$omp atomic update
grad(1, it_b) = grad(1, it_b) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_b(1)
!$omp atomic update
grad(2, it_b) = grad(2, it_b) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_b(2)
!$omp atomic update
grad(3, it_b) = grad(3, it_b) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_b(3)
end if
if (.not. sk_tc) then
!$omp atomic update
grad(1, it_c) = grad(1, it_c) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_c(1)
!$omp atomic update
grad(2, it_c) = grad(2, it_c) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_c(2)
!$omp atomic update
grad(3, it_c) = grad(3, it_c) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_c(3)
end if
if (.not. sk_td) then
!$omp atomic update
grad(1, it_d) = grad(1, it_d) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_d(1)
!$omp atomic update
grad(2, it_d) = grad(2, it_d) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_d(2)
!$omp atomic update
grad(3, it_d) = grad(3, it_d) + bds%strtork(6 + k, i) * dr3 * gt(k) * Jt_d(3)
end if
if (.not. sk_3a) then
!$omp atomic update
grad(1, ib3_a) = grad(1, ib3_a) + bds%strtork(6 + k, i) * dihef(k) * Jb3_a(1)
!$omp atomic update
grad(2, ib3_a) = grad(2, ib3_a) + bds%strtork(6 + k, i) * dihef(k) * Jb3_a(2)
!$omp atomic update
grad(3, ib3_a) = grad(3, ib3_a) + bds%strtork(6 + k, i) * dihef(k) * Jb3_a(3)
end if
if (.not. sk_3b) then
!$omp atomic update
grad(1, ib3_b) = grad(1, ib3_b) + bds%strtork(6 + k, i) * dihef(k) * Jb3_b(1)
!$omp atomic update
grad(2, ib3_b) = grad(2, ib3_b) + bds%strtork(6 + k, i) * dihef(k) * Jb3_b(2)
!$omp atomic update
grad(3, ib3_b) = grad(3, ib3_b) + bds%strtork(6 + k, i) * dihef(k) * Jb3_b(3)
end if
end do
end do
end subroutine strtor_geomgrad