Skip to content

Commit

Permalink
Fixes #1054 by ensuring thread savety of hessian calculations with GF…
Browse files Browse the repository at this point in the history 
…N-FF (#1061)

The distance array neigh%distances was not allocated thread save and was removed from the neigh class.
Now the distances are calculated from the atom coordinates and the translation vectors.
  • Loading branch information
@Thomas3R
Thomas3R authored Jul 2, 2024
1 parent 2d55ea9 commit edcfbbe
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 106 deletions.
47 changes: 27 additions & 20 deletions src/gfnff/gfnff_eg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,7 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
real(wp),allocatable :: xtmp(:)
type(tb_timer) :: timer
real(wp) :: dispthr, cnthr, repthr, hbthr1, hbthr2
real(wp) :: dist(n,n)
real(wp) :: ds(3,3)
real(wp) :: convF !convergence factor alpha, aka ewald parameter
logical, allocatable :: considered_ABH(:,:,:)
Expand Down Expand Up @@ -225,9 +226,6 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &

! get Distances between atoms for repulsion
call neigh%getTransVec(env,mol,sqrt(repthr))
if(allocated(neigh%distances)) deallocate(neigh%distances)
call getDistances(neigh%distances, mol, neigh)


g = 0
exb = 0
Expand Down Expand Up @@ -291,6 +289,18 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
sqrab(ij + i) = 0.0d0
srab(ij + i) = 0.0d0
enddo
if (mol%npbc.ne.0) then
dist = 0.0_wp
!$omp parallel do collapse(2) default(none) shared(dist,mol) &
!$omp private(i,j)
do i=1, mol%n
do j=1, mol%n
dist(j,i) = NORM2(mol%xyz(:,j)-mol%xyz(:,i))
enddo
enddo
!$omp end parallel do
endif

if (pr) call timer%measure(1)

!----------!
Expand Down Expand Up @@ -434,7 +444,7 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
if (pr) call timer%measure(4)
else
if (pr) call timer%measure(4,'EEQ energy and q')
call goed_pbc_gfnff(env,mol,accuracy.gt.1,n,at,neigh%distances(:,:,1), &
call goed_pbc_gfnff(env,mol,accuracy.gt.1,n,at,dist, &
& dfloat(ichrg),eeqtmp,cn,nlist%q,ees,solvation,param,topo, gTrans, &
& rTrans, xtmp, convF) ! without dq/dr
ees = ees*mcf_ees
Expand Down Expand Up @@ -602,7 +612,7 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
dx=xa-xyz(1,jat)-neigh%transVec(1,iTr)
dy=ya-xyz(2,jat)-neigh%transVec(2,iTr)
dz=za-xyz(3,jat)-neigh%transvec(3,iTr)
rab=neigh%distances(jat,iat,iTr)
rab=NORM2(xyz(:,iat)-(xyz(:,jat)+neigh%transVec(:,iTr)))
r2=rab**2
ati=at(iat)
atj=at(jat)
Expand Down Expand Up @@ -734,8 +744,6 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
call neigh%getTransVec(env,mol,sqrt(hbthr2))
if (pr) call timer%measure(10,'HB/XB (incl list setup)')
if (update.or.require_update) then
if(allocated(neigh%distances)) deallocate(neigh%distances)
call getDistances(neigh%distances, mol, neigh)
call gfnff_hbset(n,at,xyz,topo,neigh,nlist,hbthr1,hbthr2)
end if

Expand Down Expand Up @@ -781,12 +789,11 @@ subroutine gfnff_eg(env,mol,pr,n,ichrg,at,xyz,sigma,g,etot,res_gff, &
nbk=0
iTr=0 ! nbk is the first neighbor of k !
atnb=0
call neigh%jth_nb(nbk,1,k,iTr)
call neigh%jth_nb(n,xyz,nbk,1,k,iTr)
! get iTrC and cycle if out of cutoff (the cutoff used in last getTransVec call)
iTrDum=neigh%fTrSum(iTr,iTrk)
if(iTrDum.eq.-1.or.iTrDum.gt.neigh%nTrans) cycle
! get number neighbors of neighbor of k !
if(nbk.le.0 .or. nbk.gt.n) write(*,*) 'nbk=',nbk
if(nbk.ne.0) then
nbnbk=sum(neigh%nb(neigh%numnb,nbk,:))
atnb=at(nbk)
Expand Down Expand Up @@ -2283,7 +2290,7 @@ subroutine abhgfnff_eg2new(n,A,B,H,iTrA,iTrB,nbb,at,xyz,q,sqrab, &
! Neighbours of B
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when inb is shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when inb is shifted to iTr
! compute distances
vecDum = neigh%transVec(:,iTr)+neigh%transVec(:,iTrB)
dranb(1:3,i) = (xyz(1:3,A)+neigh%transVec(1:3,iTrA)) -(xyz(1:3,inb)+vecDum)
Expand Down Expand Up @@ -2456,7 +2463,7 @@ subroutine abhgfnff_eg2new(n,A,B,H,iTrA,iTrB,nbb,at,xyz,q,sqrab, &
gdr(1:3,H) = gdr(1:3,H) + gh(1:3)
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
gdr(1:3,inb) = gdr(1:3,inb) + gnb(1:3,i)
end do

Expand All @@ -2466,7 +2473,7 @@ subroutine abhgfnff_eg2new(n,A,B,H,iTrA,iTrB,nbb,at,xyz,q,sqrab, &
sigma=sigma+mcf_ehb*spread(gh,1,3)*spread(xyz(:,H),2,3)
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
vecDum = neigh%transVec(:,iTr)+neigh%transVec(:,iTrB)
sigma=sigma+mcf_ehb*spread(gnb(:,i),1,3)*spread(xyz(:,inb)+vecDum,2,3)
enddo
Expand Down Expand Up @@ -2543,7 +2550,7 @@ subroutine abhgfnff_eg2_rnr(n,A,B,H,iTrA,iTrB,at,xyz,q,sqrab,srab,energy,gdr,par
! Neighbours of B
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
! compute distances
vTrinb=neigh%transVec(:,iTr)+neigh%transVec(:,iTrB)
dranb(1:3,i) = (xyz(1:3,A)+neigh%transVec(1:3,iTrA)) -(xyz(1:3,inb)+vTrinb)
Expand Down Expand Up @@ -2753,7 +2760,7 @@ subroutine abhgfnff_eg2_rnr(n,A,B,H,iTrA,iTrB,at,xyz,q,sqrab,srab,energy,gdr,par
gdr(1:3,H) = gdr(1:3,H) + gh(1:3)
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
gdr(1:3,inb) = gdr(1:3,inb) + gnb(1:3,i) - gnb_lp(1:3)/dble(nbb)
end do

Expand All @@ -2764,7 +2771,7 @@ subroutine abhgfnff_eg2_rnr(n,A,B,H,iTrA,iTrB,at,xyz,q,sqrab,srab,energy,gdr,par
sigma=sigma+mcf_ehb*spread(gnb_lp,1,3)*spread(xyz(:,B)+neigh%transVec(1:3,iTrB),2,3)
do i=1,nbb
inb=0; iTr=0 ! jth_nb output
call neigh%jth_nb(inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
call neigh%jth_nb(n,xyz,inb,i,B,iTr) ! inb is the i-th nb of B when shifted to iTr
vTrinb=neigh%transVec(:,iTr)+neigh%transVec(:,iTrB)
sigma=sigma+mcf_ehb*spread(gnb(:,i),1,3)*spread(xyz(:,inb)+vTrinb,2,3)
sigma=sigma-mcf_ehb*spread(gnb_lp(1:3)/dble(nbb),1,3)*spread(xyz(:,inb)+vTrinb,2,3)
Expand Down Expand Up @@ -3267,8 +3274,8 @@ subroutine batmgfnff_eg(n,iat,jat,kat,iTrj,iTrk,at,xyz,q,sqrab,srab,energy,g,ds,
real*8 r2ij,r2jk,r2ik,c9,mijk,imjk,ijmk,rijk3,ang,angr9,rav3
real*8 rij(3),rik(3),rjk(3),ri(3),rj(3),rk(3),drij,drik,drjk,dang,ff,fi,fj,fk,fqq
parameter (fqq=3.0d0)
integer linij,linik,linjk,lina,i,j,iTrDum,dm1,dm2
lina(i,j)=min(i,j)+max(i,j)*(max(i,j)-1)/2
integer :: linij,linik,linjk,lina,i,j,iTrDum,dm1,dm2
lina(i,j)=min(i,j)+max(i,j)*(max(i,j)-1)/2

fi=(1.d0-fqq*q(iat))
fi=min(max(fi,-4.0d0),4.0d0)
Expand All @@ -3278,13 +3285,13 @@ subroutine batmgfnff_eg(n,iat,jat,kat,iTrj,iTrk,at,xyz,q,sqrab,srab,energy,g,ds,
fk=min(max(fk,-4.0d0),4.0d0)
ff=fi*fj*fk ! charge term
c9=ff*param%zb3atm(at(iat))*param%zb3atm(at(jat))*param%zb3atm(at(kat)) ! strength of interaction
r2ij=neigh%distances(jat,iat,iTrj)**2
r2ik=neigh%distances(kat,iat,iTrk)**2
r2ij=NORM2(xyz(:,iat)-(xyz(:,jat)+neigh%transVec(:,iTrj)))**2
r2ik=NORM2(xyz(:,iat)-(xyz(:,kat)+neigh%transVec(:,iTrk)))**2
iTrDum=neigh%fTrSum(neigh%iTrNeg(iTrj),iTrk)
if(iTrDum.le.0.or.iTrDum.gt.neigh%numctr) then
r2jk=NORM2((xyz(:,kat)+neigh%transVec(:,iTrk))-(xyz(:,jat)+neigh%transVec(:,iTrj)))**2
else
r2jk=neigh%distances(kat,jat,iTrDum)**2 !use adjusted iTr since jat is actually shifted
r2jk=NORM2(xyz(:,jat)-(xyz(:,kat)+neigh%transVec(:,iTrDum)))**2
endif
mijk=-r2ij+r2jk+r2ik
imjk= r2ij-r2jk+r2ik
Expand Down
38 changes: 17 additions & 21 deletions src/gfnff/gfnff_ini.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,19 +206,15 @@ integer function itabrow6(i)
sqrab = 0
!Get all distances
call neigh%getTransVec(env,mol,sqrt(hbthr2))
if (.not. allocated(neigh%distances)) then
call getDistances(neigh%distances, mol, neigh)
endif

! Transfer calculated distances in central cell to rab and sqrab
! want to change to neigh%distances
do i=1, mol%n
ati=mol%at(i)
kk=i*(i-1)/2
do j=1, i-1
atj=mol%at(j)
k=kk+j
rab(k) = neigh%distances(j,i,1)
rab(k) = NORM2(mol%xyz(:,i)-mol%xyz(:,j))
sqrab(k) = rab(k)**2
if(rab(k).lt.1.d-3) then
write(env%unit,*) i,j,ati,atj,rab(k)
Expand Down Expand Up @@ -404,7 +400,7 @@ integer function itabrow6(i)
nn =sum(neigh%nb(neigh%numnb,i,:))
if(nn.eq.0) cycle
ip =piadr2(i)
call neigh%jth_nb(ji,1,i,iTrtmp) ! ji is the first nb of i in cell iTr
call neigh%jth_nb(mol%n,mol%xyz,ji,1,i,iTrtmp) ! ji is the first nb of i in cell iTr
nh =0
nm =0
do iTr=1, neigh%numctr
Expand Down Expand Up @@ -858,7 +854,7 @@ integer function itabrow6(i)
end do
do i = 1,mol%n
! get first nb, iTr expected to be irrelevant since same in each cell
call neigh%jth_nb(nn,1,i,iTr) ! R - N/C/O - H ! terminal H
call neigh%jth_nb(mol%n,mol%xyz,nn,1,i,iTr) ! R - N/C/O - H ! terminal H
if(nn.le.0) cycle ! cycle if there is no nb
topo%hbaci(i)=param%xhaci(mol%at(i)) ! only first neighbor of interest
if (amideH(mol%n,mol%at,topo%hyb,neigh%numnb,neigh%numctr,neigh%nb,piadr2,i,neigh)) &
Expand All @@ -876,7 +872,7 @@ integer function itabrow6(i)
if(mol%at(i).ne.1) cycle
if(topo%hyb(i).eq.1) cycle ! exclude bridging hydrogens from HB correction
ff=gen%hqabthr
call neigh%jth_nb(j,1,i,iTr) ! get j, first neighbor of H when j is shifted to iTr
call neigh%jth_nb(mol%n,mol%xyz,j,1,i,iTr) ! get j, first neighbor of H when j is shifted to iTr
if(j.le.0) cycle
if(mol%at(j).gt.10) ff=ff-0.20 ! H on heavy atoms may be negatively charged
if(mol%at(j).eq.6.and.topo%hyb(j).eq.3) ff=ff+0.05 ! H on sp3 C must be really positive 0.05
Expand Down Expand Up @@ -1038,7 +1034,7 @@ integer function itabrow6(i)
ja=piadr4(jj)
if(ia.gt.0.and.ja.gt.0) then
!dum=1.d-9*rab(lin(ii,jj)) ! distort so that Huckel for e.g. COT localizes to right bonds
dum=1.d-9*neigh%distances(jj,ii,iTr) ! distort so that Huckel for e.g. COT localizes to right bonds
dum=1.d-9*NORM2(mol%xyz(:,ii)-(mol%xyz(:,jj)+neigh%transVec(:,iTr))) ! distort so that Huckel for e.g. COT localizes to right bonds
dum=sqrt(gen%hoffdiag(mol%at(ii))*gen%hoffdiag(mol%at(jj)))-dum ! better than arithmetic
dum2=gen%hiter
if(topo%hyb(ii).eq.1) dum2=dum2*gen%htriple ! triple bond is different
Expand Down Expand Up @@ -2065,16 +2061,16 @@ integer function itabrow6(i)
endif
deallocate(locarr)
! sort atoms according to distance to central atom such that the same inversion angle def. always results
sdum3(1)=neigh%distances(jj,i,iTrj)
sdum3(2)=neigh%distances(kk,i,iTrk)
sdum3(3)=neigh%distances(ll,i,iTrl)
sdum3(1)=NORM2(mol%xyz(:,i)-(mol%xyz(:,jj)+neigh%transVec(:,iTrj)))
sdum3(2)=NORM2(mol%xyz(:,i)-(mol%xyz(:,kk)+neigh%transVec(:,iTrk)))
sdum3(3)=NORM2(mol%xyz(:,i)-(mol%xyz(:,ll)+neigh%transVec(:,iTrl)))
ind3(1)=jj
ind3(2)=kk
ind3(3)=ll
call ssort(3,sdum3,ind3)
sdum3(1)=neigh%distances(jj,i,iTrj) ! need old indices for sorting iTr's
sdum3(2)=neigh%distances(kk,i,iTrk)
sdum3(3)=neigh%distances(ll,i,iTrl)
sdum3(1)=NORM2(mol%xyz(:,i)-(mol%xyz(:,jj)+neigh%transVec(:,iTrj)))
sdum3(2)=NORM2(mol%xyz(:,i)-(mol%xyz(:,kk)+neigh%transVec(:,iTrk)))
sdum3(3)=NORM2(mol%xyz(:,i)-(mol%xyz(:,ll)+neigh%transVec(:,iTrl)))
jj=ind3(1) ! assign sorted indices
kk=ind3(2)
ll=ind3(3)
Expand Down Expand Up @@ -2163,7 +2159,7 @@ integer function itabrow6(i)
if (mol%at(i).eq.6.and.sum(neigh%nb(neigh%numnb,i,:)).eq.2) then
do j=1, 2
nbi=0
call neigh%jth_nb(nbi,j,i,iTr) ! nbi is the jth nb of i in cell iTr
call neigh%jth_nb(mol%n,mol%xyz,nbi,j,i,iTr) ! nbi is the jth nb of i in cell iTr
if (nbi.eq.0) cycle
if (mol%at(nbi).eq.6.and.sum(neigh%nb(neigh%numnb,nbi,:)).eq.2) then
nn = nn + 1
Expand Down Expand Up @@ -2222,7 +2218,7 @@ subroutine specialTorsList(nst, mol, topo, neigh, sTorsList)
! carbon with two neighbors bonded to other carbon* with two neighbors
if (mol%at(i).eq.6.and.sum(neigh%nb(neigh%numnb,i,:)).eq.2) then
do j=1, 2
call neigh%jth_nb(nbi,j,i,iTr) ! nbi is the jth nb of i in cell iTr
call neigh%jth_nb(mol%n,mol%xyz,nbi,j,i,iTr) ! nbi is the jth nb of i in cell iTr
if (nbi.eq.0.or.iTr.eq.0) cycle
if (mol%at(nbi).eq.6.and.sum(neigh%nb(neigh%numnb,nbi,:)).eq.2) then ! *other carbon
! check carbon triple bond distance
Expand All @@ -2231,14 +2227,14 @@ subroutine specialTorsList(nst, mol, topo, neigh, sTorsList)
! check C2 and C3
do k=1, 2 ! C2 is other nb of Ci
nbk=0
call neigh%jth_nb(nbk,k,i,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
call neigh%jth_nb(mol%n,mol%xyz,nbk,k,i,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
if (nbk.ne.nbi.and.nbk.ne.0.and.mol%at(nbk).eq.6) then
jj=nbk ! C2 index
endif
enddo
do k=1, 2 ! C3 is other nb of Cnbi
nbk=0
call neigh%jth_nb(nbk,k,nbi,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
call neigh%jth_nb(mol%n,mol%xyz,nbk,k,nbi,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
if (nbk.ne.i.and.nbk.ne.0.and.mol%at(nbk).eq.6) then
kk=nbk ! C3 index
endif
Expand All @@ -2255,7 +2251,7 @@ subroutine specialTorsList(nst, mol, topo, neigh, sTorsList)
! on atom sorting in input file !!! The last one in file.
do k=1, sum(neigh%nb(neigh%numnb,jj,:))
nbk=0
call neigh%jth_nb(nbk,k,jj,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
call neigh%jth_nb(mol%n,mol%xyz,nbk,k,jj,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
if (topo%hyb(k).eq.2.and.mol%at(k).eq.6.and.sum(neigh%nb(neigh%numnb,k,:)).eq.3.and. &
& nbk.ne.i.and.nbk.ne.0) then
ii=nbk
Expand All @@ -2266,7 +2262,7 @@ subroutine specialTorsList(nst, mol, topo, neigh, sTorsList)
! on atom sorting in input file !!! The last one in file.
do k=1, sum(neigh%nb(neigh%numnb,kk,:))
nbk=0
call neigh%jth_nb(nbk,k,jj,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
call neigh%jth_nb(mol%n,mol%xyz,nbk,k,jj,iTr) ! nbk is the jth nb of i in cell iTr .ne.nbi
if (topo%hyb(k).eq.2.and.mol%at(k).eq.6.and.sum(neigh%nb(neigh%numnb,i,:)).eq.3.and. &
& nbk.ne.nbi.and.nbk.ne.0) then
ll=nbk
Expand Down
Loading

0 comments on commit edcfbbe

Please sign in to comment.