compute forces from long-range atoms on quantum fragment and vice versa
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(atom_type), | intent(in), | dimension(:) | :: | atoms | list of atoms in the long-range part |
|
| type(quantum_fragment_type), | intent(in) | :: | quantum_fragment | the quantum fragment |
||
| integer, | intent(in) | :: | atom_start | index of the first atom that current process will handle |
||
| integer, | intent(in) | :: | atom_end | index of the last atom that the current process will handle |
subroutine compute_lr_forces(atoms, quantum_fragment, &
atom_start, atom_end)
!> list of atoms in the long-range part
type(atom_type), dimension(:), intent(in) :: atoms
!> the quantum fragment
type(quantum_fragment_type), intent(in) :: quantum_fragment
!> index of the first atom that current process will handle
integer, intent(in) :: atom_start
!> index of the last atom that the current process will handle
integer, intent(in) :: atom_end
integer :: i, j
integer :: x, y, z
integer :: tx, ty, tz
integer :: mx, my, mz
integer :: rx, ry, rz
integer :: fx, fy, fz
integer :: mi, fi, ti
integer :: t_start, t_end, m_start, m_end
integer :: rank_f, rank_t, rank_m
integer :: order, expansion_order
real(dp) :: coeff, temp
real(dp), dimension(3) :: r
real(dp), dimension(3) :: force
real(dp), dimension(3) :: expansion_origin
real(dp), dimension(:), allocatable :: multipoles
real(dp), dimension(:), allocatable :: tensor, tensor_sum
real(dp), dimension(:), allocatable :: polytensor, polytensor_sum
call timer_start("compute_lr_forces")
if (atom_start > atom_end) then
call handle_error(SEVERITY_FATAL, &
TYPE_INCORRECT_ARG, &
"atom_start should be smaller than atom_end", &
__FILE__, __LINE__)
endif
expansion_origin = quantum_fragment%nuclear_multipoles%origin
expansion_order = quantum_fragment%nuclear_multipoles%order
allocate(multipoles(polytensor_size(expansion_order)))
multipoles = quantum_fragment%nuclear_multipoles%values + &
quantum_fragment%electronic_multipoles%values
allocate(tensor(tensor_size(expansion_order + 1)))
allocate(tensor_sum(tensor_size(expansion_order + 1)))
allocate(polytensor(polytensor_size(expansion_order)))
allocate(polytensor_sum(polytensor_size(expansion_order)))
!$OMP PARALLEL DO private(i, r, force, order, m_start, m_end, &
!$OMP& rank_t, rank_m, rank_f, fx, fy, fz, fi, tx, ty, tz, ti, &
!$OMP& mx, my, mz, mi, coeff, tensor)
do i = atom_start, atom_end
r = expansion_origin - atoms(i)%coordinate
force = 0.0_dp
m_end = 0
do order = 0, expansion_order
m_start = m_end + 1
m_end = m_start + tensor_size(order) - 1
call interaction_tensor(r, tensor(1:tensor_size(order+1)), order, 1)
rank_t = tensor_rank(tensor(1:tensor_size(order+1)))
rank_m = tensor_rank(multipoles(m_start:m_end))
rank_f = tensor_rank(force)
do fx = rank_f, 0, -1
do fy = rank_f - fx, 0, -1
fz = rank_f - fx - fy
fi = convert_index(fx, fy, fz)
do tx = rank_t, fx, -1
do ty = rank_t - tx, fy, -1
tz = rank_t - tx - ty
if (tz < fz) cycle
ti = convert_index(tx, ty, tz)
mx = tx - fx
my = ty - fy
mz = tz - fz
mi = convert_index(mx, my, mz)
coeff = real((-1)**(tx+ty+tz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
force(fi) = force(fi) + coeff * tensor(ti) * multipoles(m_start-1+mi)
end do
end do
end do
end do
end do
atoms(i)%force = atoms(i)%force - atoms(i)%charge * force
end do
!$OMP END PARALLEL DO
tensor_sum = 0.0_dp
!$OMP PARALLEL DO private(i, r, tensor) reduction(+:tensor_sum)
do i = atom_start, atom_end
r = expansion_origin - atoms(i)%coordinate
call interaction_tensor(r, tensor, expansion_order, 1)
tensor_sum = tensor_sum + atoms(i)%charge * tensor
end do
!$OMP END PARALLEL DO
force = 0.0_dp
m_start = polytensor_size(expansion_order - 1) + 1
m_end = m_start + tensor_size(expansion_order) - 1
rank_t = tensor_rank(tensor_sum)
rank_m = tensor_rank(multipoles(m_start:m_end))
rank_f = tensor_rank(force)
do fx = rank_f, 0, -1
do fy = rank_f - fx, 0, -1
fz = rank_f - fx - fy
fi = convert_index(fx, fy, fz)
do tx = rank_t, fx, -1
do ty = rank_t - tx, fy, -1
tz = rank_t - tx - ty
if (tz < fz) cycle
ti = convert_index(tx, ty, tz)
mx = tx - fx
my = ty - fy
mz = tz - fz
mi = convert_index(mx, my, mz)
coeff = real((-1)**(mx+my+mz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
force(fi) = force(fi) + coeff * tensor_sum(ti) * multipoles(m_start-1+mi)
end do
end do
end do
end do
force = force / real(quantum_fragment%num_nuclei, dp)
do j = 1, quantum_fragment%num_nuclei
quantum_fragment%nuclei(j)%force = quantum_fragment%nuclei(j)%force - force
end do
polytensor = 0.0_dp
polytensor_sum = 0.0_dp
!$OMP PARALLEL DO private(i, r, order, t_start, t_end, polytensor) &
!$OMP& reduction(+:polytensor_sum)
do i = atom_start, atom_end
r = expansion_origin - atoms(i)%coordinate
t_end = 1
do order = 1, expansion_order
t_start = t_end + 1
t_end = t_start + tensor_size(order) - 1
call interaction_tensor(r, polytensor(t_start:t_end), order, 0)
end do
polytensor_sum = polytensor_sum + atoms(i)%charge * polytensor
end do
!$OMP END PARALLEL DO
!$OMP PARALLEL DO private(j, r, multipoles, order, rx, ry, rz, temp, &
!$OMP& force, t_start, t_end, m_start, m_end, rank_t, rank_m, rank_f, &
!$OMP& fx, fy, fz, fi, tx, ty, tz, ti, mx, my, mz, mi, coeff)
do j = 1, quantum_fragment%num_nuclei
r = quantum_fragment%nuclei(j)%coordinate - expansion_origin
multipoles = 0.0_dp
mi = 1
do order = 0, expansion_order - 1
do rx = order, 0, -1
do ry = order - rx, 0, -1
rz = order - rx - ry
temp = quantum_fragment%nuclei(j)%charge
do x = 1, rx
temp = temp * r(1)
end do
do y = 1, ry
temp = temp * r(2)
end do
do z = 1, rz
temp = temp * r(3)
end do
multipoles(mi) = multipoles(mi) + temp
mi = mi + 1
end do
end do
end do
force = 0.0_dp
t_end = 1
m_end = 0
do order = 1, expansion_order
t_start = t_end + 1
t_end = t_start + tensor_size(order) - 1
m_start = m_end + 1
m_end = m_start + tensor_size(order-1) - 1
rank_t = tensor_rank(polytensor_sum(t_start:t_end))
rank_m = tensor_rank(multipoles(m_start:m_end))
rank_f = tensor_rank(force)
do fx = rank_f, 0, -1
do fy = rank_f - fx, 0, -1
fz = rank_f - fx - fy
fi = convert_index(fx, fy, fz)
do tx = rank_t, fx, -1
do ty = rank_t - tx, fy, -1
tz = rank_t - tx - ty
if (tz < fz) cycle
ti = convert_index(tx, ty, tz)
mx = tx - fx
my = ty - fy
mz = tz - fz
mi = convert_index(mx, my, mz)
coeff = real((-1)**(tx+ty+tz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
force(fi) = force(fi) + coeff * polytensor_sum(t_start-1+ti) * multipoles(m_start-1+mi)
end do
end do
end do
end do
end do
quantum_fragment%nuclei(j)%force = quantum_fragment%nuclei(j)%force - force
end do
!$OMP END PARALLEL DO
deallocate(multipoles)
deallocate(tensor)
deallocate(tensor_sum)
deallocate(polytensor)
deallocate(polytensor_sum)
call timer_stop
end subroutine compute_lr_forces