Skip to content
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

Jump to bottom

Fixes #1054 by ensuring thread savety of hessian calculations with GFN-FF #1061

Merged
merged 2 commits into from
Jul 2, 2024
Merged

Fixes #1054 by ensuring thread savety of hessian calculations with GFN-FF #1061

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9ff3d5a
Fix thread savety of hessian calculations with GFN-FF
Thomas3R Jul 2, 2024
55070b3
Update gfnff_ini2.f90
Thomas3R Jul 2, 2024
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
47 changes: 27 additions & 20 deletions src/gfnff/gfnff_eg.f90
View file
Open in desktop
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
Open in desktop
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
Loading