| 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 angtor_geomgrad(bds, grad)
use mod_jacobian_mat, only: simple_angle_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), da1, da2, angle1, angle2, f1, f2, f3, &
Jt_a(3), Jt_b(3), Jt_c(3), Jt_d(3), &
Ja1_a(3), Ja1_b(3), Ja1_c(3), &
Ja2_a(3), Ja2_b(3), Ja2_c(3)
integer(ip) :: i, j, k, ia1, ia2, &
it_a, it_b, it_c, it_d, &
ia1_a, ia1_b, ia1_c, &
ia2_a, ia2_b, ia2_c
logical :: sk_ta, sk_tb, sk_tc, sk_td, &
sk_1a, sk_1b, sk_1c, &
sk_2a, sk_2b, sk_2c
if(.not. bds%use_angtor) return
!$omp parallel do default(shared) &
!$omp private(thet, gt, dihef, da1, da2, angle1, angle2, f1, f2, f3, Jt_a, Jt_b) &
!$omp private(Jt_c, Jt_d, Ja1_a, Ja1_b, Ja1_c) &
!$omp private(Ja2_a, Ja2_b, Ja2_c, i, j, k, ia1, ia2) &
!$omp private(it_a, it_b, it_c, it_d, ia1_a, ia1_b, ia1_c, ia2_a, ia2_b, ia2_c) &
!$omp private(sk_ta, sk_tb, sk_tc, sk_td, sk_1a, sk_1b, sk_1c, sk_2a, sk_2b, sk_2c)
do i=1, bds%nangtor
! Atoms that define the dihedral angle
it_a = bds%angtorat(1,i)
it_b = bds%angtorat(2,i)
it_c = bds%angtorat(3,i)
it_d = bds%angtorat(4,i)
ia1 = bds%angtor_a(1,i)
ia1_a = bds%angleat(1,ia1)
ia1_b = bds%angleat(2,ia1)
ia1_c = bds%angleat(3,ia1)
ia2 = bds%angtor_a(2,i)
ia2_a = bds%angleat(1,ia2)
ia2_b = bds%angleat(2,ia2)
ia2_c = bds%angleat(3,ia2)
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(ia1_a)
sk_1b = bds%top%frozen(ia1_b)
sk_1c = bds%top%frozen(ia1_c)
sk_2a = bds%top%frozen(ia2_a)
sk_2b = bds%top%frozen(ia2_b)
sk_2c = bds%top%frozen(ia2_c)
if(sk_ta .and. sk_tb .and. sk_tc .and. sk_td .and. &
sk_1a .and. sk_1b .and. sk_1c .and. &
sk_2a .and. sk_2b .and. sk_2c) cycle
else
sk_ta = .false.
sk_tb = .false.
sk_tc = .false.
sk_td = .false.
sk_1a = .false.
sk_1b = .false.
sk_1c = .false.
sk_2a = .false.
sk_2b = .false.
sk_2c = .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 simple_angle_jacobian(bds%top%cmm(:,ia1_a), &
bds%top%cmm(:,ia1_b), &
bds%top%cmm(:,ia1_c), &
angle1, Ja1_a, Ja1_b, Ja1_c)
call simple_angle_jacobian(bds%top%cmm(:,ia2_a), &
bds%top%cmm(:,ia2_b), &
bds%top%cmm(:,ia2_c), &
angle2, Ja2_a, Ja2_b, Ja2_c)
da1 = angle1 - bds%eqangle(ia1)
da2 = angle2 - bds%eqangle(ia2)
do k = 1, 3
if(.not.(sk_ta .and. sk_tb .and. sk_tc .and. sk_td)) &
f1 = (bds%angtork(k, i) * da1 + bds%angtork(3+k,i) * da2) * gt(k)
if(.not.(sk_1a .and. sk_1b .and. sk_1c)) &
f2 = bds%angtork(k, i) * dihef(k)
if(.not.(sk_2a .and. sk_2b .and. sk_2c)) &
f3 = bds%angtork(3+k, i) * dihef(k)
if (.not. sk_ta) then
!$omp atomic update
grad(1, it_a) = grad(1, it_a) + f1 * Jt_a(1)
!$omp atomic update
grad(2, it_a) = grad(2, it_a) + f1 * Jt_a(2)
!$omp atomic update
grad(3, it_a) = grad(3, it_a) + f1 * Jt_a(3)
end if
if (.not. sk_tb) then
!$omp atomic update
grad(1, it_b) = grad(1, it_b) + f1 * Jt_b(1)
!$omp atomic update
grad(2, it_b) = grad(2, it_b) + f1 * Jt_b(2)
!$omp atomic update
grad(3, it_b) = grad(3, it_b) + f1 * Jt_b(3)
end if
if (.not. sk_tc) then
!$omp atomic update
grad(1, it_c) = grad(1, it_c) + f1 * Jt_c(1)
!$omp atomic update
grad(2, it_c) = grad(2, it_c) + f1 * Jt_c(2)
!$omp atomic update
grad(3, it_c) = grad(3, it_c) + f1 * Jt_c(3)
end if
if (.not. sk_td) then
!$omp atomic update
grad(1, it_d) = grad(1, it_d) + f1 * Jt_d(1)
!$omp atomic update
grad(2, it_d) = grad(2, it_d) + f1 * Jt_d(2)
!$omp atomic update
grad(3, it_d) = grad(3, it_d) + f1 * Jt_d(3)
end if
if (.not. sk_1a) then
!$omp atomic update
grad(1, ia1_a) = grad(1, ia1_a) + f2 * Ja1_a(1)
!$omp atomic update
grad(2, ia1_a) = grad(2, ia1_a) + f2 * Ja1_a(2)
!$omp atomic update
grad(3, ia1_a) = grad(3, ia1_a) + f2 * Ja1_a(3)
end if
if (.not. sk_1b) then
!$omp atomic update
grad(1, ia1_b) = grad(1, ia1_b) + f2 * Ja1_b(1)
!$omp atomic update
grad(2, ia1_b) = grad(2, ia1_b) + f2 * Ja1_b(2)
!$omp atomic update
grad(3, ia1_b) = grad(3, ia1_b) + f2 * Ja1_b(3)
end if
if (.not. sk_1c) then
!$omp atomic update
grad(1, ia1_c) = grad(1, ia1_c) + f2 * Ja1_c(1)
!$omp atomic update
grad(2, ia1_c) = grad(2, ia1_c) + f2 * Ja1_c(2)
!$omp atomic update
grad(3, ia1_c) = grad(3, ia1_c) + f2 * Ja1_c(3)
end if
if (.not. sk_2a) then
!$omp atomic update
grad(1, ia2_a) = grad(1, ia2_a) + f3 * Ja2_a(1)
!$omp atomic update
grad(2, ia2_a) = grad(2, ia2_a) + f3 * Ja2_a(2)
!$omp atomic update
grad(3, ia2_a) = grad(3, ia2_a) + f3 * Ja2_a(3)
end if
if (.not. sk_2b) then
!$omp atomic update
grad(1, ia2_b) = grad(1, ia2_b) + f3 * Ja2_b(1)
!$omp atomic update
grad(2, ia2_b) = grad(2, ia2_b) + f3 * Ja2_b(2)
!$omp atomic update
grad(3, ia2_b) = grad(3, ia2_b) + f3 * Ja2_b(3)
end if
if (.not. sk_2c) then
!$omp atomic update
grad(1, ia2_c) = grad(1, ia2_c) + f3 * Ja2_c(1)
!$omp atomic update
grad(2, ia2_c) = grad(2, ia2_c) + f3 * Ja2_c(2)
!$omp atomic update
grad(3, ia2_c) = grad(3, ia2_c) + f3 * Ja2_c(3)
end if
end do
end do
end subroutine angtor_geomgrad