Compute angle-bending terms of the potential energy function.
Simple angle terms are computed according to the formula:
Out-of plane angle are more complex. First, central atom has to be
a trigonal center, the other two atoms together with the auxliary
atom (that is the remaining one connected to the trigonal center)
define the projection plane. During the first run the auxiliary atom
is found and saved.
Then, the trigonal center is projected on the plane defined by the
other three atoms, and the angle is the one defined by the projection
(which is the vertex, and the other two atoms -- the auxiliary is
excluded). Then the same formula used for simple angle terms is used.
Trigonal center
Normal vector of the projection plane
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_bonded_type), | intent(in) | :: | bds | |||
| real(kind=rp), | intent(inout) | :: | V |
Bond potential, result will be added to V |
subroutine angle_potential(bds, V)
!! Compute angle-bending terms of the potential energy function.
!! Simple angle terms are computed according to the formula:
!! \[U_{angle} = \sum_i k_i \Delta \theta_i^2 \large(1 +
!! \sum_{j=1}^4 k^{(j+2)} \Delta \theta_i^j \large)\]
!! \[\Delta \theta_i = \theta_i - \theta^{(eq)}_i\]
!! Out-of plane angle are more complex. First, central atom has to be
!! a trigonal center, the other two atoms together with the auxliary
!! atom (that is the remaining one connected to the trigonal center)
!! define the projection plane. During the first run the auxiliary atom
!! is found and saved.
!! Then, the trigonal center is projected on the plane defined by the
!! other three atoms, and the angle is the one defined by the projection
!! (which is the vertex, and the other two atoms -- the auxiliary is
!! excluded). Then the same formula used for simple angle terms is used.
use mod_constants, only: eps_rp
implicit none
type(ommp_bonded_type), intent(in) :: bds
! Bonded potential data structure
real(rp), intent(inout) :: V
!! Bond potential, result will be added to V
integer(ip) :: i
real(rp) :: l1, l2, dr1(3), dr2(3), thet, d_theta
real(rp), dimension(3) :: v_dist, plv1, plv2, pln, a, b, c, prj_b, aux
if(.not. bds%use_angle) return
!$omp parallel do default(shared) schedule(static) reduction(+:V) &
!$omp private(i,dr1,dr2,l1,l2,thet,d_theta,a,b,c,aux,plv1,plv2,pln,v_dist,prj_b)
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
dr1 = bds%top%cmm(:, bds%angleat(1,i)) - bds%top%cmm(:, bds%angleat(2,i))
dr2 = bds%top%cmm(:, bds%angleat(3,i)) - bds%top%cmm(:, bds%angleat(2,i))
l1 = sqrt(dot_product(dr1, dr1))
l2 = sqrt(dot_product(dr2, dr2))
thet = acos(dot_product(dr1, dr2)/(l1*l2))
d_theta = thet-bds%eqangle(i)
V = V + bds%kangle(i) * d_theta**2 * (1.0 + bds%angle_cubic*d_theta &
+ bds%angle_quartic*d_theta**2 + bds%angle_pentic*d_theta**3 &
+ bds%angle_sextic*d_theta**4)
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
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))
plv1 = a - aux
plv2 = c - aux
pln(1) = plv1(2)*plv2(3) - plv1(3)*plv2(2)
pln(2) = plv1(3)*plv2(1) - plv1(1)*plv2(3)
pln(3) = plv1(1)*plv2(2) - plv1(2)*plv2(1)
!! Normal vector of the projection plane
pln = pln / sqrt(dot_product(pln, pln))
v_dist = b - aux
prj_b = b - dot_product(v_dist, pln) * pln
dr1 = bds%top%cmm(:, bds%angleat(1,i)) - prj_b
dr2 = bds%top%cmm(:, bds%angleat(3,i)) - prj_b
l1 = sqrt(dot_product(dr1, dr1))
l2 = sqrt(dot_product(dr2, dr2))
thet = acos(dot_product(dr1, dr2)/(l1*l2))
d_theta = thet-bds%eqangle(i)
V = V + bds%kangle(i) * d_theta**2 * (1.0 + bds%angle_cubic*d_theta &
+ bds%angle_quartic*d_theta**2 + bds%angle_pentic*d_theta**3 &
+ bds%angle_sextic*d_theta**4)
end if
end do
end subroutine angle_potential