Compute the dispersion repulsion energy for the whole system using a double loop algorithm
All neighbors done!
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_nonbonded_type), | intent(in), | target | :: | vdw |
Nonbonded data structure |
|
| real(kind=rp), | intent(inout) | :: | V |
Potential, result will be added |
subroutine vdw_potential(vdw, V)
!! Compute the dispersion repulsion energy for the whole system
!! using a double loop algorithm
use mod_io, only : fatal_error
use mod_profiling, only: time_push, time_pull
use mod_constants, only: eps_rp
use mod_memory, only: mallocate, mfree
use mod_neighbor_list, only: get_ith_nl
implicit none
type(ommp_nonbonded_type), intent(in), target :: vdw
!! Nonbonded data structure
real(rp), intent(inout) :: V
!! Potential, result will be added
integer(ip) :: i, j, jc, l, ipair, ineigh, nthreads, ithread, nn
real(rp) :: eij, rij0, rij, ci(3), cj(3), s, vtmp
type(ommp_topology_type), pointer :: top
procedure(vdw_term), pointer :: vdw_func
integer(ip), allocatable :: nl_neigh(:,:)
real(rp), allocatable :: nl_r(:,:)
integer :: omp_get_num_threads, omp_get_thread_num
call time_push()
!$omp parallel
nthreads = omp_get_num_threads()
!$omp end parallel
top => vdw%top
select case(vdw%vdwtype)
case(OMMP_VDWTYPE_LJ)
vdw_func => vdw_lennard_jones
case(OMMP_VDWTYPE_BUF714)
vdw_func => vdw_buffered_7_14
case default
vdw_func => vdw_buffered_7_14
call fatal_error("Unexpected error in vdw_potential")
end select
if(vdw%use_nl) then
call mallocate('vdw_potential [rneigh]', top%mm_atoms, nthreads, nl_r)
call mallocate('vdw_potential [nl_neigh]', top%mm_atoms, nthreads, nl_neigh)
end if
!$omp parallel do default(shared) reduction(+:v) schedule(dynamic) &
!$omp private(i,j,jc,ineigh,ithread,nn,s,ci,cj,ipair,l,Eij,Rij0,Rij,vtmp)
do i=1, top%mm_atoms
ithread = omp_get_thread_num() + 1
if(abs(vdw%vdw_f(i) - 1.0_rp) < eps_rp) then
ci = top%cmm(:,i)
else
! Scale factors are used only for monovalent atoms, in that
! case the vdw center is displaced along the axis connecting
! the atom to its neighbour
if(top%conn(1)%ri(i+1) - top%conn(1)%ri(i) /= 1) then
call fatal_error("Scale factors are only expected for &
&monovalent atoms")
end if
ineigh = top%conn(1)%ci(top%conn(1)%ri(i))
ci = top%cmm(:,ineigh) + (top%cmm(:,i) - top%cmm(:,ineigh)) &
* vdw%vdw_f(i)
endif
! If neighbor list are enabled get the one for the current
if(vdw%use_nl) call get_ith_nl(vdw%nl, i, top%cmm, &
nl_neigh(:,ithread), &
nl_r(:,ithread), nn)
do jc=1, top%mm_atoms
! If the two atoms aren't neighbors, just skip the loop
if(vdw%use_nl) then
if(jc > nn) exit !! All neighbors done!
j = nl_neigh(jc,ithread)
if(j <= i) cycle
else
! Skip all iteration with j <= i
if(jc > i) then
j = jc
else
cycle
end if
end if
! Compute the screening factor for this pair
s = 1.0_rp
do ineigh=1,4
! Look if j is at distance ineigh from i
if(any(top%conn(ineigh)%ci(top%conn(ineigh)%ri(i): &
top%conn(ineigh)%ri(i+1)-1) == j)) then
s = vdw%vdw_screening(ineigh)
! Exit the loop
exit
end if
end do
if(s > eps_rp) then
ipair = -1
do l=1, vdw%npair
if((vdw%vdw_pair_mask_a(i,l) .and. vdw%vdw_pair_mask_b(j,l)) .or. &
(vdw%vdw_pair_mask_a(j,l) .and. vdw%vdw_pair_mask_b(i,l))) then
ipair = l
exit
end if
end do
if(ipair > 0) then
Rij0 = vdw%vdw_pair_r(ipair)
eij = vdw%vdw_pair_e(ipair)
else
Rij0 = get_Rij0(vdw, i, j)
eij = get_eij(vdw, i, j)
end if
if(abs(vdw%vdw_f(j) - 1.0_rp) < eps_rp) then
cj = top%cmm(:,j)
else
! Scale factors are used only for monovalent atoms, in that
! case the vdw center is displaced along the axis connecting
! the atom to its neighbour
if(top%conn(1)%ri(j+1) - top%conn(1)%ri(j) /= 1) then
call fatal_error("Scale factors are only expected &
& for monovalent atoms")
end if
ineigh = top%conn(1)%ci(top%conn(1)%ri(j))
cj = top%cmm(:,ineigh) + &
(top%cmm(:,j) - top%cmm(:,ineigh)) * vdw%vdw_f(j)
endif
Rij = norm2(ci-cj)
if(Rij < eps_rp) then
call fatal_error("Requesting non-bonded potential for two atoms &
&placed in the same point, this could be &
&an internal bug or a problem in your input &
&file, please check.")
end if
vtmp = 0.0_rp
call vdw_func(Rij, Rij0, Eij, vtmp)
v = v + vtmp*s
end if
end do
end do
if(vdw%use_nl) then
call mfree('vdw_potential [rneigh]', nl_r)
call mfree('vdw_potential [nl_neigh]', nl_neigh)
end if
call time_pull('VdW potential calculation')
end subroutine vdw_potential