| 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 bond_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) :: ca(3), cb(3), J_a(3), J_b(3), l, dl, g
use_cubic = (abs(bds%bond_cubic) > eps_rp)
use_quartic = (abs(bds%bond_quartic) > eps_rp)
if(.not. bds%use_bond) return
if(.not. use_cubic .and. .not. use_quartic) then
! This is just a regular harmonic potential
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,ia,ib,sk_a,sk_b,ca,cb,dl,l,g,J_a,J_b)
do i=1, bds%nbond
ia = bds%bondat(1,i)
ib = bds%bondat(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
ca = bds%top%cmm(:,ia)
cb = bds%top%cmm(:,ib)
call Rij_jacobian(ca, cb, l, J_a, J_b)
dl = l - bds%l0bond(i)
g = 2 * bds%kbond(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) schedule(dynamic) &
!$omp private(i,ia,ib,sk_a,sk_b,ca,cb,dl,l,g,J_a,J_b)
do i=1, bds%nbond
ia = bds%bondat(1,i)
ib = bds%bondat(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
ca = bds%top%cmm(:,ia)
cb = bds%top%cmm(:,ib)
call Rij_jacobian(ca, cb, l, J_a, J_b)
dl = l - bds%l0bond(i)
g = 2 * bds%kbond(i) * dl * (1.0_rp + 3.0/2.0*bds%bond_cubic*dl &
+ 2.0*bds%bond_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 bond_geomgrad