Skip to content

Commit

Permalink
Remove the slower calculation from vectorized Coulomb real
Browse files Browse the repository at this point in the history
  • Loading branch information
@mewall @cnegre
mewall authored and cnegre committed Apr 11, 2024
1 parent d09d418 commit d4cd4bd
Showing 1 changed file with 0 additions and 177 deletions.
177 changes: 0 additions & 177 deletions src/latte_mods/coulomb_latte_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -880,183 +880,6 @@ subroutine get_ewald_list_real_dcalc_vect(spindex,splist,coordinates,charges,hub
enddo
!$omp end parallel do


!coul_pot_r = 0.0_dp
!coul_forces_r = 0.0_dp

!$omp parallel do default(none) private(i) &
!$omp private(fcoul,coulombv,dcoulombv,dforce) &
!$omp private(ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse) &
!$omp private(tj,tj2,tj3,tj4,tj6,ti2mtj2,sa,sb,sc,sd,se,sf) &
!$omp private(ra,rb,nni,dr,rab,magr,magr2,j) &
!$omp private(dc,z,numrep_erfc,ca,force,expti,exptj,tj2mti2,rmod) &
!$omp shared(nats,hubbardu,spindex,coordinates,sqrtpi,keconst,Lx,Ly,Lz ) &
!$omp shared(nrnnlist,coulcut,nnType,tfact,nnIx,nnIy,nnIz,splist) &
!$omp shared(sa_mat,sb_mat,sc_mat,sd_mat,se_mat,sf_mat) &
!$omp shared(e_mat,e1_mat,e2_mat,f_mat,f1_mat,f2_mat) &
!$omp shared(coul_forces_r, coul_pot_r, calpha, charges, calpha2)
do i =1,nats

fcoul = 0.0_dp
coulombv = 0.0_dp

ti = tfact*hubbardu(spindex(i));

ti2 = ti*ti;
ti3 = ti2*ti;
ti4 = ti2*ti2;
ti6 = ti4*ti2;

ssa = ti;
ssb = ti3/48.0_dp;
ssc = 3.0_dp*ti2/16.0_dp;
ssd = 11.0_dp*ti/16.0_dp;
sse = 1.0_dp;

ra = coordinates(:,i);

do nni = 1,nrnnlist(i)

j = nnType(nni,i);

if(allocated(nnIx))then
Rb(1) = coordinates(1,j) + nnIx(nni,i)*Lx
Rb(2) = coordinates(2,j) + nnIy(nni,i)*Ly
Rb(3) = coordinates(3,j) + nnIz(nni,i)*Lz
rab = rb-ra
rmod = prg_norm2(rab)

else

Rb(1) = coordinates(1,j)
Rb(2) = coordinates(2,j)
Rb(3) = coordinates(3,j)

rab = rb-ra

rmod = prg_norm2(rab)

if(rmod > coulcut)then
rab(1) = modulo((Rb(1) - Ra(1) + Lx/2.0_dp),Lx) - Lx/2.0_dp
rab(2) = modulo((Rb(2) - Ra(2) + Ly/2.0_dp),Ly) - Ly/2.0_dp
rab(3) = modulo((Rb(3) - Ra(3) + Lz/2.0_dp),Lz) - Lz/2.0_dp
endif
endif

dr = norm2(rab)

magr = dr
magr2 = dr*dr

if (dr <= coulcut .and. dr > 1e-12) then
tj = tfact*hubbardu(spindex(j))
dc = rab/dr
z = abs(calpha*magr)
numrep_erfc = erfc(z)
ca = numrep_erfc/magr
coulombv = coulombv + charges(j)*ca
ca = ca + 2.0_dp*calpha*exp( -calpha2*magr2 )/sqrtpi
force = -keconst*charges(i)*charges(j)*ca/magr
expti = exp(-ti*magr)

if (splist(spindex(i)) == splist(spindex(j)))then
dcoulombv = e1_mat(nni,i)
!dcoulombv = - charges(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
coulombv = coulombv + dcoulombv
dforce = f1_mat(nni,i)
!dforce = (keconst*charges(i)*charges(j)*expti)*((sse/magr2 - 2*ssb*magr - ssc) +&
! ssa*(ssb*magr2 + ssc*magr + ssd + sse/magr))
force = force + dforce
! if(abs(e1_mat(nni,i)-dcoulombv).gt.1.0D-12)then
! write(*,*)"e1_mat != dcoulombv",nni,i,e1_mat(nni,i),dcoulombv
! endif
! if(abs(f1_mat(nni,i)-dforce).gt.1.0D-12)then
! write(*,*)"f1_mat != dforce",nni,i,f1_mat(nni,i),dforce
! endif
else
tj2 = tj*tj
tj3 = tj2*tj
tj4 = tj2*tj2
tj6 = tj4*tj2
exptj = exp( -tj*magr )
ti2mtj2 = ti2 - tj2
tj2mti2 = -ti2mtj2
sa = ti
sb = tj4*ti/(2.0_dp * ti2mtj2 * ti2mtj2)
sc = (tj6 - 3.0_dp*tj4*ti2)/(ti2mtj2 * ti2mtj2 * ti2mtj2)
sd = tj
se = ti4*tj/(2.0_dp * tj2mti2 * tj2mti2)
sf = (ti6 - 3.0_dp*ti4*tj2)/(tj2mti2 * tj2mti2 * tj2mti2)
! if(sa_mat(nni,i).ne.sa)then
! write(*,*)"sa_mat != sa",nni,i,sa_mat(nni,i),sa
! else
! write(*,*)"sa_mat == sa",nni,i
! endif
! if(sb_mat(nni,i).ne.sb)then
! write(*,*)"sb_mat != sb",nni,i,sb_mat(nni,i),sb
! else
! write(*,*)"sb_mat == sb",nni,i
! endif
! if(sc_mat(nni,i).ne.sc)then
! write(*,*)"sc_mat != sc",nni,i,sc_mat(nni,i),sc
! else
! write(*,*)"sc_mat == sc",nni,i
! endif
! if(sd_mat(nni,i).ne.sd)then
! write(*,*)"sd_mat != sd",nni,i,sd_mat(nni,i),sd
! else
! write(*,*)"sd_mat == sd",nni,i
! endif
! if(se_mat(nni,i).ne.se)then
! write(*,*)"se_mat != se",nni,i,se_mat(nni,i),se
! else
! write(*,*)"se_mat == se",nni,i
! endif
! if(sf_mat(nni,i).ne.sf)then
! write(*,*)"sf_mat != sf",nni,i,sf_mat(nni,i),sf
! else
! write(*,*)"sf_mat == sf",nni,i
! endif
dcoulombv = e2_mat(nni,i)
!dcoulombv = - (charges(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
coulombv = coulombv + dcoulombv
dforce = f2_mat(nni,i)
!dforce = keconst*charges(i)*charges(j)*((expti*(sa*(sb - (sc/magr)) - (sc/magr2))) +&
! (exptj*(sd*(se - (sf/magr)) - (sf/magr2))))
force = force + dforce
! if(abs(e2_mat(nni,i)-dcoulombv).gt.1.0D-12)then
! write(*,*)"e2_mat != dcoulombv",nni,i,e2_mat(nni,i),dcoulombv
! endif
! if(abs(f2_mat(nni,i)-dforce).gt.1.0D-12)then
! write(*,*)"f2_mat != dforce",nni,i,f2_mat(nni,i),dforce
! endif
endif

!if(abs(dforce-f_mat(nnI,i)).gt.1.0D-12)then
! write(*,*)"Forces differ for ",nnI,i,f_mat(nnI,i),f1_mat(nnI,i),f2_mat(nnI,i),force
! stop
!endif
fcoul = fcoul + dc*force

endif

enddo
!$omp critical
!if(any(abs(fcoul-coul_forces_r(:,i)).gt.1.0D-12))then
! write(*,*)"Forces differ for atom index ",i,fcoul(1),coul_forces_r(1,i)
!endif
if(abs(keconst*coulombv-coul_pot_r(i)).gt.1.0D-12)then
write(*,*)"Potentials differ for atom index ",i,keconst*coulombv,coul_pot_r(i)
endif
!coul_forces_r(:,i) = fcoul
!coul_pot_r(i) = coulombv
!$omp end critical
enddo
!$omp end parallel do


!coul_pot_r = keconst*coul_pot_r

end subroutine get_ewald_list_real_dcalc_vect

!> This routine computes the reciprocal space contribution of the Ewald summation.
Expand Down

0 comments on commit d4cd4bd

Please sign in to comment.