Compute torsion potential
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 imptorsion_geomgrad(bds, grad)
!! Compute torsion potential
use mod_jacobian_mat, only: 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, g, J_a(3), J_b(3), J_c(3), J_d(3)
integer(ip) :: i, j, ia, ib, ic, id
logical :: sk_a, sk_b, sk_c, sk_d
if (.not. bds%use_imptorsion) return
!$omp parallel do default(shared) &
!$omp private(i, ia, ib, ic, id, sk_a, sk_b, sk_c, sk_d, j, thet, J_a, J_b, J_c, J_d, g)
do i = 1, bds%nimptorsion
! Atoms that define the dihedral angle
ia = bds%imptorsionat(1, i)
ib = bds%imptorsionat(2, i)
ic = bds%imptorsionat(3, i)
id = bds%imptorsionat(4, 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)
sk_d = bds%top%frozen(id)
if (sk_a .and. sk_b .and. sk_c .and. sk_d) cycle
else
sk_a = .false.
sk_b = .false.
sk_c = .false.
sk_d = .false.
end if
call torsion_angle_jacobian(bds%top%cmm(:, ia), &
bds%top%cmm(:, ib), &
bds%top%cmm(:, ic), &
bds%top%cmm(:, id), &
thet, J_a, J_b, J_c, J_d)
do j = 1, 3
if (bds%imptorsn(j, i) < 1) exit
g = -real(bds%imptorsn(j, i)) * sin(real(bds%imptorsn(j, i)) * thet &
- bds%imptorsphase(j, i)) &
* bds%imptorsamp(j, i)
if (.not. sk_a) then
!$omp atomic update
grad(1, ia) = grad(1, ia) + J_a(1) * g
!$omp atomic update
grad(2, ia) = grad(2, ia) + J_a(2) * g
!$omp atomic update
grad(3, ia) = grad(3, ia) + J_a(3) * g
end if
if (.not. sk_b) then
!$omp atomic update
grad(1, ib) = grad(1, ib) + J_b(1) * g
!$omp atomic update
grad(2, ib) = grad(2, ib) + J_b(2) * g
!$omp atomic update
grad(3, ib) = grad(3, ib) + J_b(3) * g
end if
if (.not. sk_c) then
!$omp atomic update
grad(1, ic) = grad(1, ic) + J_c(1) * g
!$omp atomic update
grad(2, ic) = grad(2, ic) + J_c(2) * g
!$omp atomic update
grad(3, ic) = grad(3, ic) + J_c(3) * g
end if
if (.not. sk_d) then
!$omp atomic update
grad(1, id) = grad(1, id) + J_d(1) * g
!$omp atomic update
grad(2, id) = grad(2, id) + J_d(2) * g
!$omp atomic update
grad(3, id) = grad(3, id) + J_d(3) * g
end if
end do
end do
end subroutine imptorsion_geomgrad