Skip to content

Implement parallelization for hydrogen bond list setup in gfnff_hbset #1240

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

Implement parallelization for hydrogen bond list setup in gfnff_hbset #1240

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
7b3626f
Add OpenMP parallelization to HB setup in gfnff_hbset
Thomas3R Mar 27, 2025
6ec47c6
Merge branch 'devel' of github.com:Thomas3R/xtb into main
Thomas3R Mar 27, 2025
310a0c1
Optimize distance calculation and readability in gfnff_hbset
Thomas3R Mar 28, 2025
3ec282f
Update if conditions in gfnff_hbset
Thomas3R Mar 28, 2025
fef5d59
Use mctc_dot in gfnff_hbset
Thomas3R Mar 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
148 changes: 87 additions & 61 deletions src/gfnff/gfnff_ini2.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -763,6 +763,7 @@ end function chkrng

subroutine gfnff_hbset(n,at,xyz,topo,neigh,nlist,hbthr1,hbthr2)
use xtb_mctc_accuracy, only : wp
use xtb_mctc_blas_level1, only : mctc_dot
use xtb_gfnff_param
implicit none
type(TGFFTopology), intent(in) :: topo
Expand All @@ -774,85 +775,110 @@ subroutine gfnff_hbset(n,at,xyz,topo,neigh,nlist,hbthr1,hbthr2)
real(wp), intent(in) :: hbthr1, hbthr2

integer i,j,k,nh,ia,ix,lin,ij,inh,jnh, iTri,iTrj,iTrDum
real(wp) rab,rmsd, rih,rjh
real(wp) rab,rmsd, rih,rjh, vec_ab(3), vec_ih(3), vec_jh(3)
logical ijnonbond,free
! local copies for omp parallelization
integer :: neigh_nTrans, neigh_numctr, nlist_nhb1, nlist_nhb2

neigh_nTrans = neigh%nTrans
neigh_numctr = neigh%numctr

rmsd = sqrt(sum((xyz-nlist%hbrefgeo)**2))/dble(n)

if(rmsd.lt.1.d-6 .or. rmsd.gt. 0.3d0) then ! update list if first call or substantial move occured
nlist%nhb1=0
nlist%nhb2=0
nlist_nhb1=0
nlist_nhb2=0
! loop over hb-relevant AB atoms
!$omp parallel do default(none) private(i, j, k, nh, ix, iTri, iTrj, iTrDum, rab, &
!$omp rih, rjh,ijnonbond,free, vec_ab, vec_ih, vec_jh ) &
!$omp firstprivate (nlist_nhb1, nlist_nhb2, neigh_numctr, neigh_nTrans, hbthr1, hbthr2) &
!$omp shared(nlist, topo, neigh, xyz)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe schedule(dynamic) improves performance here. The loop work could be different due to the conditions in the nested loops.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will test the influence of the schedule in the implementation in gfnff_hbset0. The parallelization in gfnff_hbset has little impact in comparison to gfnff_hbset0, so it is easier for me to check there.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better to use collapse(3) with evaluation of i/j variables in inner loops to increase OpenMP iteration space. Playing with second arg of dynamic scheduling may also help: I usually use 32 or 64.

do ix=1,topo%nathbAB
i=topo%hbatABl(1,ix)
j=topo%hbatABl(2,ix)
do iTri=1, neigh%nTrans ! go through i shifts
do iTrj=1, neigh%nTrans ! go through j shifts
! get adjustet iTr -> for use of neigh% distances and bpair with two shifted atoms
iTrDum=neigh%fTrSum(neigh%iTrNeg(iTri),iTrj)
if(iTrDum.gt.neigh%nTrans.or.iTrDum.lt.-1.or.iTrDum.eq.0) cycle ! cycle nonsense
rab=NORM2((xyz(:,i)+neigh%transVec(:,iTri))-(xyz(:,j)+neigh%transVec(:,iTrj)))**2
if(rab.gt.hbthr1) cycle
! check if ij bonded
if(iTrDum.le.neigh%numctr.and.iTrDum.gt.0) then
ijnonbond=neigh%bpair(j,i,iTrDum).ne.1
else
! i and j are not in neighboring cells
ijnonbond=.true.
endif
! loop over relevant H atoms
do k=1,topo%nathbH
free=.true. ! tripplet not assigned yet
nh =topo%hbatHl(1,k) ! nh always in central cell
! distances for non-cov bonded case
rih=NORM2(xyz(:,nh)-(xyz(:,i)+neigh%transVec(:,iTri)))**2
rjh=NORM2(xyz(:,nh)-(xyz(:,j)+neigh%transVec(:,iTrj)))**2
! check if i is the bonded A
if(iTri.le.neigh%numctr) then ! nh is not shifted so bpair works without adjustment
if(neigh%bpair(i,nh,iTri).eq.1.and.ijnonbond) then
nlist%nhb2=nlist%nhb2+1
nlist%hblist2(1,nlist%nhb2)=i
nlist%hblist2(2,nlist%nhb2)=j
nlist%hblist2(3,nlist%nhb2)=nh
nlist%hblist2(4,nlist%nhb2)=iTri
nlist%hblist2(5,nlist%nhb2)=iTrj
free=.false. ! not available for nhb1 !!!
endif
endif
! check if j is the bonded A
if(iTrj.le.neigh%numctr.and.free) then
if(neigh%bpair(j,nh,iTrj).eq.1.and.ijnonbond) then
nlist%nhb2=nlist%nhb2+1
nlist%hblist2(1,nlist%nhb2)=j
nlist%hblist2(2,nlist%nhb2)=i
nlist%hblist2(3,nlist%nhb2)=nh
nlist%hblist2(4,nlist%nhb2)=iTrj
nlist%hblist2(5,nlist%nhb2)=iTri
free=.false. ! not available for nhb1 !!!
endif
endif
! check for non-cov bonded A
if(rab+rih+rjh.lt.hbthr2.and.free) then ! sum of rAB,rAH,rBH is below threshold
nlist%nhb1=nlist%nhb1+1
nlist%hblist1(1,nlist%nhb1)=i
nlist%hblist1(2,nlist%nhb1)=j
nlist%hblist1(3,nlist%nhb1)=nh
nlist%hblist1(4,nlist%nhb1)=iTri
nlist%hblist1(5,nlist%nhb1)=iTrj
endif
enddo ! k: relevant H atoms
enddo ! iTrj
enddo ! iTri
i=topo%hbatABl(1,ix)
j=topo%hbatABl(2,ix)
do iTri=1, neigh_nTrans ! go through i shifts
do iTrj=1, neigh_nTrans ! go through j shifts
Comment on lines +801 to +802
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice naming of iterators. :|

! get adjustet iTr -> for use of neigh% distances and bpair with two shifted atoms
iTrDum=neigh%fTrSum(neigh%iTrNeg(iTri),iTrj)
! iTrDum=-1 is valid here because we are only interested if ij is a bond (ijnonbond)
! However, iTrDum is not and should not be used as an index without excluding -1 value
if(iTrDum.gt.neigh%nTrans.or.iTrDum.lt.-1.or.iTrDum.eq.0) cycle
vec_ab = (xyz(:,i)+neigh%transVec(:,iTri))-(xyz(:,j)+neigh%transVec(:,iTrj))
rab = mctc_dot(vec_ab, vec_ab) ! square of distance AB
if(rab.gt.hbthr1) cycle
! check if ij bonded
if(iTrDum.le.neigh_numctr.and.iTrDum.gt.0) then
ijnonbond=neigh%bpair(j,i,iTrDum).ne.1
else
! i and j are not in neighboring cells
ijnonbond=.true.
endif
! loop over relevant H atoms
do k=1,topo%nathbH
free=.true. ! tripplet not assigned yet
nh =topo%hbatHl(1,k) ! nh always in central cell
! distances for non-cov bonded case
vec_ih = xyz(:,nh)-(xyz(:,i)+neigh%transVec(:,iTri))
rih = mctc_dot(vec_ih, vec_ih) ! square of distance iH
vec_jh = xyz(:,nh)-(xyz(:,j)+neigh%transVec(:,iTrj))
rjh = mctc_dot(vec_jh, vec_jh) ! square of distance jH
Comment on lines +824 to +826
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is better to use dot_product intrinsic (which can be optimized by compiler) rather than using external procedure mctc_dot

! check if i is the bonded A
if(iTri.le.neigh_numctr) then ! nh is not shifted so bpair works without adjustment
if(neigh%bpair(i,nh,iTri).eq.1.and.ijnonbond) then
!$omp atomic capture
nlist%nhb2=nlist%nhb2+1
nlist_nhb2 = nlist%nhb2
!$omp end atomic
nlist%hblist2(1,nlist_nhb2)=i
nlist%hblist2(2,nlist_nhb2)=j
nlist%hblist2(3,nlist_nhb2)=nh
nlist%hblist2(4,nlist_nhb2)=iTri
nlist%hblist2(5,nlist_nhb2)=iTrj
free=.false. ! not available for nhb1 !!!
endif
endif
! check if j is the bonded A
if(iTrj.le.neigh_numctr.and.free) then
if(neigh%bpair(j,nh,iTrj).eq.1.and.ijnonbond) then
!$omp atomic capture
nlist%nhb2=nlist%nhb2+1
nlist_nhb2 = nlist%nhb2
!$omp end atomic
nlist%hblist2(1,nlist_nhb2)=j
nlist%hblist2(2,nlist_nhb2)=i
nlist%hblist2(3,nlist_nhb2)=nh
nlist%hblist2(4,nlist_nhb2)=iTrj
nlist%hblist2(5,nlist_nhb2)=iTri
Comment on lines +849 to +853
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here you may have data race with nlist_nhb2 variable. On x86/x86_64 it works fine, but AArch64 can say "Hi" for you :)

free=.false. ! not available for nhb1 !!!
endif
endif
! check for non-cov bonded A
if(rab+rih+rjh.lt.hbthr2.and.free) then ! sum of rAB,rAH,rBH is below threshold
!$omp atomic capture
nlist%nhb1=nlist%nhb1+1
nlist_nhb1 = nlist%nhb1
!$omp end atomic
nlist%hblist1(1,nlist_nhb1)=i
nlist%hblist1(2,nlist_nhb1)=j
nlist%hblist1(3,nlist_nhb1)=nh
nlist%hblist1(4,nlist_nhb1)=iTri
nlist%hblist1(5,nlist_nhb1)=iTrj
endif
enddo ! k: relevant H atoms
enddo ! iTrj
enddo ! iTri
enddo ! ix: relevant AB atoms
!$omp end parallel do

! for nxb list only i is not shifted
nlist%nxb =0
do ix=1,topo%natxbAB
i =topo%xbatABl(1,ix) ! A
j =topo%xbatABl(2,ix) ! B
iTrj=topo%xbatABl(4,ix) ! iTrB
if(iTrj.gt.neigh%nTrans.or.iTrj.lt.-1.or.iTrj.eq.0) cycle ! cycle nonsense
if(iTrj.gt.neigh_nTrans.or.iTrj.lt.-1.or.iTrj.eq.0) cycle ! cycle nonsense
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why was it changed?

rab=NORM2(xyz(:,j)-xyz(:,i)+neigh%transVec(:,iTrj))**2
if(rab.gt.hbthr2)cycle
nlist%nxb=nlist%nxb+1
Expand Down