| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_bonded_type), | intent(in) | :: | bds |
Bonded potential data structure |
||
| real(kind=rp), | intent(inout) | :: | grad(3,bds%top%mm_atoms) |
Gradients of bond stretching terms of potential energy |
subroutine urey_geomgrad(bds, grad)
use mod_constants, only : eps_rp
use mod_jacobian_mat, only: Rij_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 :: i, ia, ib
logical(lp) :: use_cubic, use_quartic
logical :: sk_a, sk_b
real(rp) :: l, dl, J_a(3), J_b(3), g
if(.not. bds%use_urey) return
use_cubic = (abs(bds%urey_cubic) > eps_rp)
use_quartic = (abs(bds%urey_quartic) > eps_rp)
if(.not. use_cubic .and. .not. use_quartic) then
! This is just a regular harmonic potential
!$omp parallel do default(shared) &
!$omp private(i,ia,ib,sk_a,sk_b,l,dl,g,J_a,J_b)
do i=1, bds%nurey
ia = bds%ureyat(1,i)
ib = bds%ureyat(2,i)
if(bds%top%use_frozen) then
sk_a = bds%top%frozen(ia)
sk_b = bds%top%frozen(ib)
if(sk_a .and. sk_b) cycle
else
sk_a = .false.
sk_b = .false.
end if
call Rij_jacobian(bds%top%cmm(:,ia), &
bds%top%cmm(:,ib), &
l, J_a, J_b)
dl = l - bds%l0urey(i)
g = 2 * bds%kurey(i) * dl
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
end do
else
!$omp parallel do default(shared) &
!$omp private(i,ia,ib,sk_a,sk_b,l,dl,g,J_a,J_b)
do i=1, bds%nurey
ia = bds%ureyat(1,i)
ib = bds%ureyat(2,i)
if(bds%top%use_frozen) then
sk_a = bds%top%frozen(ia)
sk_b = bds%top%frozen(ib)
if(sk_a .and. sk_b) cycle
else
sk_a = .false.
sk_b = .false.
end if
call Rij_jacobian(bds%top%cmm(:,ia), &
bds%top%cmm(:,ib), &
l, J_a, J_b)
dl = l - bds%l0urey(i)
g = 2 * bds%kurey(i) * dl * (1.0 &
+ 3.0/2.0 * bds%urey_cubic*dl &
+ 2.0 * bds%urey_quartic*dl**2)
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
end do
end if
end subroutine urey_geomgrad