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 tortor_geomgrad(bds, grad)
!! Compute torsion potential
use mod_utils, only: compute_bicubic_interp
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) :: thetx, thety, vtt, dvttdx, dvttdy
real(rp), dimension(3) :: J1_a, J1_b, J2_b, J1_c, &
J2_c, J1_d, J2_d, J2_e
integer(ip) :: i, j, iprm, ibeg, iend, ia, ib, ic, id, ie
logical :: sk_a, sk_b, sk_c, sk_d, sk_e
if(.not. bds%use_tortor) return
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,iprm,ibeg,j,iend,ia,ib,ic,id,ie,sk_a,sk_b,sk_c,sk_d,sk_e) &
!$omp private(thetx,thety,J1_a,J1_b,J1_c,J1_d,J2_b,J2_c,J2_d,J2_e,vtt,dvttdx,dvttdy)
do i=1, bds%ntortor
! Atoms that defines the two angles
iprm = bds%tortorprm(i)
ibeg = 1
do j=1, iprm-1
ibeg = ibeg + bds%ttmap_shape(1,j)*bds%ttmap_shape(2,j)
end do
iend = ibeg + bds%ttmap_shape(1,iprm)*bds%ttmap_shape(2,iprm) - 1
ia = bds%tortorat(1,i)
ib = bds%tortorat(2,i)
ic = bds%tortorat(3,i)
id = bds%tortorat(4,i)
ie = bds%tortorat(5,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)
sk_e = bds%top%frozen(ie)
if(sk_a .and. sk_b .and. sk_c .and. sk_d .and. sk_e) cycle
else
sk_a = .false.
sk_b = .false.
sk_c = .false.
sk_d = .false.
sk_e = .false.
end if
call torsion_angle_jacobian(bds%top%cmm(:,ia), &
bds%top%cmm(:,ib), &
bds%top%cmm(:,ic), &
bds%top%cmm(:,id), &
thetx, &
J1_a, J1_b, J1_c, J1_d)
thetx = ang_torsion(bds%top, bds%tortorat(1:4,i))
call torsion_angle_jacobian(bds%top%cmm(:,ib), &
bds%top%cmm(:,ic), &
bds%top%cmm(:,id), &
bds%top%cmm(:,ie), &
thety, &
J2_b, J2_c, J2_d, J2_e)
thety = ang_torsion(bds%top, bds%tortorat(2:5,i))
call compute_bicubic_interp(thetx, thety, vtt, &
dvttdx, dvttdy, &
bds%ttmap_shape(1,iprm), &
bds%ttmap_shape(2,iprm), &
bds%ttmap_ang1(ibeg:iend), &
bds%ttmap_ang2(ibeg:iend), &
bds%ttmap_v(ibeg:iend), &
bds%ttmap_vx(ibeg:iend), &
bds%ttmap_vy(ibeg:iend), &
bds%ttmap_vxy(ibeg:iend))
if(.not. sk_a) then
!$omp atomic update
grad(1,ia) = grad(1,ia) + J1_a(1) * dvttdx
!$omp atomic update
grad(2,ia) = grad(2,ia) + J1_a(2) * dvttdx
!$omp atomic update
grad(3,ia) = grad(3,ia) + J1_a(3) * dvttdx
end if
if(.not. sk_b) then
!$omp atomic update
grad(1,ib) = grad(1,ib) + J1_b(1) * dvttdx + J2_b(1) * dvttdy
!$omp atomic update
grad(2,ib) = grad(2,ib) + J1_b(2) * dvttdx + J2_b(2) * dvttdy
!$omp atomic update
grad(3,ib) = grad(3,ib) + J1_b(3) * dvttdx + J2_b(3) * dvttdy
end if
if(.not. sk_c) then
!$omp atomic update
grad(1,ic) = grad(1,ic) + J1_c(1) * dvttdx + J2_c(1) * dvttdy
!$omp atomic update
grad(2,ic) = grad(2,ic) + J1_c(2) * dvttdx + J2_c(2) * dvttdy
!$omp atomic update
grad(3,ic) = grad(3,ic) + J1_c(3) * dvttdx + J2_c(3) * dvttdy
end if
if(.not. sk_d) then
!$omp atomic update
grad(1,id) = grad(1,id) + J1_d(1) * dvttdx + J2_d(1) * dvttdy
!$omp atomic update
grad(2,id) = grad(2,id) + J1_d(2) * dvttdx + J2_d(2) * dvttdy
!$omp atomic update
grad(3,id) = grad(3,id) + J1_d(3) * dvttdx + J2_d(3) * dvttdy
end if
if(.not. sk_e) then
!$omp atomic update
grad(1,ie) = grad(1,ie) + J2_e(1) * dvttdy
!$omp atomic update
grad(2,ie) = grad(2,ie) + J2_e(2) * dvttdy
!$omp atomic update
grad(3,ie) = grad(3,ie) + J2_e(3) * dvttdy
end if
end do
end subroutine tortor_geomgrad