Trigonal center
| 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 angle_geomgrad(bds, grad)
use mod_jacobian_mat, only: simple_angle_jacobian, &
inplane_angle_jacobian
use mod_constants, only: eps_rp
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
real(rp) :: a(3), b(3), c(3), Ja(3), Jb(3), Jc(3), Jx(3), g, thet, &
d_theta, aux(3)
integer(ip) :: i
logical :: sk_a, sk_b, sk_c, sk_x
if(.not. bds%use_angle) return
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,sk_a,sk_b,sk_c,sk_x,a,b,c,aux,thet,d_theta,g,Ja,Jb,Jc,Jx)
do i=1, bds%nangle
if(abs(bds%kangle(i)) < eps_rp) cycle
if(bds%anglety(i) == OMMP_ANG_SIMPLE .or. &
bds%anglety(i) == OMMP_ANG_H0 .or. &
bds%anglety(i) == OMMP_ANG_H1 .or. &
bds%anglety(i) == OMMP_ANG_H2) then
if(bds%top%use_frozen) then
sk_a = bds%top%frozen(bds%angleat(1,i))
sk_b = bds%top%frozen(bds%angleat(2,i))
sk_c = bds%top%frozen(bds%angleat(3,i))
if(sk_a .and. sk_b .and. sk_c) cycle
else
sk_a = .false.
sk_b = .false.
sk_c = .false.
end if
a = bds%top%cmm(:, bds%angleat(1,i))
b = bds%top%cmm(:, bds%angleat(2,i))
c = bds%top%cmm(:, bds%angleat(3,i))
call simple_angle_jacobian(a, b, c, thet, Ja, Jb, Jc)
d_theta = thet - bds%eqangle(i)
g = bds%kangle(i) * d_theta * (2.0 &
+ 3.0 * bds%angle_cubic * d_theta &
+ 4.0 * bds%angle_quartic * d_theta**2 &
+ 5.0 * bds%angle_pentic * d_theta**3 &
+ 6.0 * bds%angle_sextic * d_theta**4)
if(.not. sk_a) then
!$omp atomic update
grad(1,bds%angleat(1,i)) = grad(1,bds%angleat(1,i)) + g * Ja(1)
!$omp atomic update
grad(2,bds%angleat(1,i)) = grad(2,bds%angleat(1,i)) + g * Ja(2)
!$omp atomic update
grad(3,bds%angleat(1,i)) = grad(3,bds%angleat(1,i)) + g * Ja(3)
end if
if(.not. sk_b) then
!$omp atomic update
grad(1,bds%angleat(2,i)) = grad(1,bds%angleat(2,i)) + g * Jb(1)
!$omp atomic update
grad(2,bds%angleat(2,i)) = grad(2,bds%angleat(2,i)) + g * Jb(2)
!$omp atomic update
grad(3,bds%angleat(2,i)) = grad(3,bds%angleat(2,i)) + g * Jb(3)
end if
if(.not. sk_c) then
!$omp atomic update
grad(1,bds%angleat(3,i)) = grad(1,bds%angleat(3,i)) + g * Jc(1)
!$omp atomic update
grad(2,bds%angleat(3,i)) = grad(2,bds%angleat(3,i)) + g * Jc(2)
!$omp atomic update
grad(3,bds%angleat(3,i)) = grad(3,bds%angleat(3,i)) + g * Jc(3)
end if
else if(bds%anglety(i) == OMMP_ANG_INPLANE .or. &
bds%anglety(i) == OMMP_ANG_INPLANE_H0 .or. &
bds%anglety(i) == OMMP_ANG_INPLANE_H1) then
if(bds%top%use_frozen) then
sk_a = bds%top%frozen(bds%angleat(1,i))
sk_b = bds%top%frozen(bds%angleat(2,i))
sk_c = bds%top%frozen(bds%angleat(3,i))
sk_x = bds%top%frozen(bds%angauxat(i))
if(sk_a .and. sk_b .and. sk_c .and. sk_x) cycle
else
sk_a = .false.
sk_b = .false.
sk_c = .false.
sk_x = .false.
end if
a = bds%top%cmm(:, bds%angleat(1,i))
b = bds%top%cmm(:, bds%angleat(2,i)) !! Trigonal center
c = bds%top%cmm(:, bds%angleat(3,i))
aux = bds%top%cmm(:, bds%angauxat(i))
call inplane_angle_jacobian(a, b, c, aux, thet, Ja, Jb, Jc, Jx)
d_theta = thet - bds%eqangle(i)
g = bds%kangle(i) * d_theta * (2.0 &
+ 3.0 * bds%angle_cubic * d_theta &
+ 4.0 * bds%angle_quartic * d_theta**2 &
+ 5.0 * bds%angle_pentic * d_theta**3 &
+ 6.0 * bds%angle_sextic * d_theta**4)
if(.not. sk_a) then
!$omp atomic update
grad(1,bds%angleat(1,i)) = grad(1,bds%angleat(1,i)) + g * Ja(1)
!$omp atomic update
grad(2,bds%angleat(1,i)) = grad(2,bds%angleat(1,i)) + g * Ja(2)
!$omp atomic update
grad(3,bds%angleat(1,i)) = grad(3,bds%angleat(1,i)) + g * Ja(3)
end if
if(.not. sk_b) then
!$omp atomic update
grad(1,bds%angleat(2,i)) = grad(1,bds%angleat(2,i)) + g * Jb(1)
!$omp atomic update
grad(2,bds%angleat(2,i)) = grad(2,bds%angleat(2,i)) + g * Jb(2)
!$omp atomic update
grad(3,bds%angleat(2,i)) = grad(3,bds%angleat(2,i)) + g * Jb(3)
end if
if(.not. sk_c) then
!$omp atomic update
grad(1,bds%angleat(3,i)) = grad(1,bds%angleat(3,i)) + g * Jc(1)
!$omp atomic update
grad(2,bds%angleat(3,i)) = grad(2,bds%angleat(3,i)) + g * Jc(2)
!$omp atomic update
grad(3,bds%angleat(3,i)) = grad(3,bds%angleat(3,i)) + g * Jc(3)
end if
if(.not. sk_x) then
!$omp atomic update
grad(1,bds%angauxat(i)) = grad(1,bds%angauxat(i)) + g * Jx(1)
!$omp atomic update
grad(2,bds%angauxat(i)) = grad(2,bds%angauxat(i)) + g * Jx(2)
!$omp atomic update
grad(3,bds%angauxat(i)) = grad(3,bds%angauxat(i)) + g * Jx(3)
end if
end if
end do
end subroutine angle_geomgrad