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

Roci #96

Closed
wants to merge 2 commits into from
Closed

Roci #96

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
76f9422
Functional changes to avoid numerical errors in epc17 low density eva…
fsmoncadaa Jun 18, 2024
e446f2a
small commit fixing a division by zero and a test that didnt pass
fsmoncadaa Nov 4, 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
34 changes: 13 additions & 21 deletions src/DFT/Functional.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -567,7 +567,7 @@ subroutine Functional_EPCEvaluate( this, mass, n, rhoE, rhoN, ec, vcE, vcN )
real(8) :: ec(*) !! Energy density - output
real(8) :: vcE(*), vcN(*) !! Potentials - output

real(8) :: a,b,c, q
real(8) :: a,b,c, prodRho
real(8) :: denominator, densityThreshold
real(8) :: v_exchange(n),va_correlation(n),vb_correlation(n)
integer :: i
Expand All @@ -588,35 +588,27 @@ subroutine Functional_EPCEvaluate( this, mass, n, rhoE, rhoN, ec, vcE, vcN )
STOP "The nuclear electron functional chosen is not implemented"
end if

!$omp parallel private(denominator)
ec(1:n)=0.0
vcE(1:n)=0.0
vcN(1:n)=0.0

! print *, "i, rhoE, rhoN, denominator, energy density, potentialE, potentialN"
!$omp parallel private(denominator,prodRho)
!$omp do schedule (dynamic)
do i = 1, n

denominator=a-b*sqrt(rhoE(i)*rhoN(i))+c*rhoE(i)*rhoN(i)

if( rhoE(i)+rhoN(i) .lt. densityThreshold ) cycle
prodRho=rhoE(i)*rhoN(i)
denominator=a-b*sqrt(prodRho)+c*prodRho
!!!Energy density
! ec(i)= -rhoE(1:n)*rhoN(1:n)/denominator(1:n)
ec(i)= -rhoN(i)/denominator

!!!Potential

if( rhoE(i)+rhoN(i) .gt. densityThreshold ) then !
vcE(i)= (rhoE(i)*rhoN(i)*(c*rhoN(i)-b*rhoN(i)/(2*sqrt(rhoE(i)*rhoN(i))))-rhoN(i)*denominator)/denominator**2
vcN(i)= (rhoN(i)*rhoE(i)*(c*rhoE(i)-b*rhoE(i)/(2*sqrt(rhoE(i)*rhoN(i))))-rhoE(i)*denominator)/denominator**2
else
vcE(i)=0.0
vcN(i)=0.0
end if

vcE(i)= rhoN(i)*(prodRho*c-sqrt(prodRho)*b/2.0-denominator)/denominator**2
vcN(i)= rhoE(i)*(prodRho*c-sqrt(prodRho)*b/2.0-denominator)/denominator**2
! write(*,"(I0.1,5ES16.6)") i, rhoE(i), rhoN(i), ec(i), vcE(i), vcN(i)
end do
!$omp end do
!$omp end parallel

! print *, "i, rhoE, rhoN, denominator, energy density, potentialE, potentialN"
! do i = 1, n
! write(*,"(I0.1,5F16.6)") i, rhoE(i), rhoN(i), ec(i), vcE(i), vcN(i)
! end do

end subroutine Functional_EPCEvaluate

subroutine Functional_IKNEvaluate( this, mass, n, rhoE, rhoN, ec, vcE, vcN )
Expand Down
22 changes: 18 additions & 4 deletions src/DFT/GridManager.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,8 +283,13 @@ subroutine GridManager_atomicOrbitals( action, type )
STOP "ERROR At DFT program, requested an unknown grid type to orbitals at GridManager"
end if

if (.not. allocated(Grid_instance(speciesID)%orbitalsWithGradient)) &
allocate(Grid_instance(speciesID)%orbitalsWithGradient(totalNumberOfContractions))
if (allocated(Grid_instance(speciesID)%orbitalsWithGradient)) then
do u=1, size(Grid_instance(speciesID)%orbitalsWithGradient(:))
call Matrix_destructor(Grid_instance(speciesID)%orbitalsWithGradient(u))
end do
deallocate(Grid_instance(speciesID)%orbitalsWithGradient)
end if
allocate(Grid_instance(speciesID)%orbitalsWithGradient(totalNumberOfContractions))

do u=1, totalNumberOfContractions
call Matrix_Constructor( Grid_instance(speciesID)%orbitalsWithGradient(u), int(gridSize,8), int(4,8), 0.0_8)
Expand Down Expand Up @@ -314,7 +319,7 @@ subroutine GridManager_atomicOrbitals( action, type )
call Matrix_constructor( auxMatrix(2), int(gridSize,8), int(numberOfCartesiansOrbitals,8), 0.0_8) !d orbital/dx
call Matrix_constructor( auxMatrix(3), int(gridSize,8), int(numberOfCartesiansOrbitals,8), 0.0_8) !d orbital/dy
call Matrix_constructor( auxMatrix(4), int(gridSize,8), int(numberOfCartesiansOrbitals,8), 0.0_8) !d orbital/dz

call ContractedGaussian_getGradientAtGrid( MolecularSystem_instance%species(speciesID)%particles(g)%basis%contraction(i), &
Grid_instance(speciesID)%points, gridSize, auxMatrix(1), auxMatrix(2), auxMatrix(3), auxMatrix(4))

Expand Down Expand Up @@ -375,8 +380,8 @@ subroutine GridManager_getDensityGradientAtGrid( speciesID, densityMatrix, densi
integer :: i, j, u, g
integer :: ii, jj, v, gg
integer :: s, ss
real(8) :: sum
integer :: numberOfContractions
real(8) :: auxreal

integer :: n, nproc
real :: time1,time2
Expand Down Expand Up @@ -451,6 +456,15 @@ subroutine GridManager_getDensityGradientAtGrid( speciesID, densityMatrix, densi
call Vector_Destructor( nodeGradientInGrid(n,3))
end do

!!Check for negative values, which should be numerical mistakes
auxreal=1.0E8
do point = 1 , gridSize
if(densityInGrid%values(point).lt.auxreal) auxreal=densityInGrid%values(point)
if(densityInGrid%values(point).lt.0.0) densityInGrid%values(point)=0.0
end do
if(-auxreal.gt.CONTROL_instance%NUCLEAR_ELECTRON_DENSITY_THRESHOLD) print *, "found significative negative density, up to", auxreal, ", for species", speciesID


time2=omp_get_wtime()
! write(*,"(A,F10.3,A4)") "**getDensityGradientAtGrid:", time2-time1 ," (s)"

Expand Down
Loading
Loading