LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Fixed floating-point exceptions

Browse Source
This commit is contained in:
Anthony Scemama 2020-05-25 18:11:27 +02:00
parent 75891f14b7
commit a90b446beb
7 changed files with 93 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
bin/qp_reset
View File

@ -112,7 +112,7 @@ qp_edit --check ${ezfio}
if [[ $mos -eq 1 ]] ; then if [[ $mos -eq 1 ]] ; then
qp set mo_two_e_ints io_mo_two_e_integrals None qp set mo_two_e_ints io_mo_two_e_integrals None
qp set mo_one_e_ints io_mo_integrals_e_n None qp set mo_one_e_ints io_mo_integrals_n_e None
qp set mo_one_e_ints io_mo_integrals_kinetic None qp set mo_one_e_ints io_mo_integrals_kinetic None
qp set mo_one_e_ints io_mo_integrals_pseudo None qp set mo_one_e_ints io_mo_integrals_pseudo None
qp set mo_one_e_ints io_mo_one_e_integrals None qp set mo_one_e_ints io_mo_one_e_integrals None

2
config/gfortran_debug.cfg
View File

@ -51,7 +51,7 @@ FCFLAGS : -Ofast
# -g : Extra debugging information # -g : Extra debugging information
# #
[DEBUG] [DEBUG]
FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
# OpenMP flags # OpenMP flags
################# #################

2
config/travis.cfg
View File

@ -53,7 +53,7 @@ FCFLAGS : -Ofast -fimplicit-none
# -g : Extra debugging information # -g : Extra debugging information
# #
[DEBUG] [DEBUG]
FCFLAGS : -Ofast -fcheck=all -g -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant FCFLAGS : -Ofast -fcheck=all -g -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
# OpenMP flags # OpenMP flags

6
src/ao_one_e_ints/pseudopot.f90
View File

@ -1907,13 +1907,17 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
sum=s_q_0 sum=s_q_0
if (q>300) then if (q>300) then
stop 'pseudopot.f90 : q > 300' stop 'pseudopot.f90 : q > 200'
endif endif
qk = dble(q) qk = dble(q)
two_qkmp1 = 2.d0*(qk+mk)+1.d0 two_qkmp1 = 2.d0*(qk+mk)+1.d0
do k=0,q-1 do k=0,q-1
s_q_k = two_qkmp1*qk*inverses(k)*s_q_k s_q_k = two_qkmp1*qk*inverses(k)*s_q_k
if (s_q_k < 1.d-32) then
s_q_k = 0.d0
exit
endif
sum=sum+s_q_k sum=sum+s_q_k
two_qkmp1 = two_qkmp1-2.d0 two_qkmp1 = two_qkmp1-2.d0
qk = qk-1.d0 qk = qk-1.d0

13
src/hartree_fock/10.hf.bats
View File

@ -5,12 +5,13 @@ source $QP_ROOT/quantum_package.rc
function run() { function run() {
thresh=1.e-7 thresh=1.e-6
test_exe scf || skip test_exe scf || skip
qp set_file $1 qp set_file $1
qp edit --check qp edit --check
qp reset --mos qp reset --mos
qp run scf qp set scf_utils n_it_scf_max 50
qp run scf
# qp set_frozen_core # qp set_frozen_core
energy="$(ezfio get hartree_fock energy)" energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh eq $energy $2 $thresh
@ -39,7 +40,7 @@ function run() {
} }
@test "SO" { # 0.539000 5.70403s @test "SO" { # 0.539000 5.70403s
run so.ezfio -25.7175263371941 run so.ezfio -25.7175270084056
} }
@test "HCO" { # 0.636700 1.55279s @test "HCO" { # 0.636700 1.55279s
@ -107,13 +108,13 @@ function run() {
} }
@test "C2H2" { # 19.599000 37.7923s @test "C2H2" { # 19.599000 37.7923s
run c2h2.ezfio -12.12144019495306 run c2h2.ezfio -12.12144044853196
} }
@test "SiH3" { # 20.316100 54.0861s @test "SiH3" { # 20.316100 54.0861s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
run sih3.ezfio -5.455398769158780 run sih3.ezfio -5.455400439077580
} }
@test "OH" { # 32.042200 1.36478m @test "OH" { # 32.042200 1.36478m
@ -130,6 +131,6 @@ function run() {
@test "SO2" { # 71.894900 3.22567m @test "SO2" { # 71.894900 3.22567m
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
run so2.ezfio -41.55800190733211 run so2.ezfio -41.55800401346361
} }

24
src/nuclei/nuclei.irp.f
View File

@ -53,18 +53,20 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ]
'================','============','============','============','============' '================','============','============','============','============'
write(6,'(A)') '' write(6,'(A)') ''
double precision :: dist_min, x, y, z if (nucl_num > 1) then
dist_min = huge(1.d0) double precision :: dist_min, x, y, z
do i=1,nucl_num dist_min = huge(1.d0)
do j=i+1,nucl_num do i=1,nucl_num
x = nucl_coord(i,1)-nucl_coord(j,1) do j=i+1,nucl_num
y = nucl_coord(i,2)-nucl_coord(j,2) x = nucl_coord(i,1)-nucl_coord(j,1)
z = nucl_coord(i,3)-nucl_coord(j,3) y = nucl_coord(i,2)-nucl_coord(j,2)
dist_min = min(x*x + y*y + z*z, dist_min) z = nucl_coord(i,3)-nucl_coord(j,3)
dist_min = min(x*x + y*y + z*z, dist_min)
enddo
enddo enddo
enddo write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', &
write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', & dsqrt(dist_min)*a0,' Angstrom'
dsqrt(dist_min)*a0,' Angstrom' endif
endif endif

120
src/scf_utils/roothaan_hall_scf.irp.f
View File

@ -196,6 +196,13 @@ END_DOC
double precision,allocatable :: scratch(:,:) double precision,allocatable :: scratch(:,:)
integer :: i,j,k,i_DIIS,j_DIIS integer :: i,j,k,i_DIIS,j_DIIS
double precision :: rcond, ferr, berr
integer, allocatable :: iwork(:)
integer :: lwork
if (dim_DIIS < 4) then
return
endif
allocate( & allocate( &
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), & B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), &
@ -239,77 +246,80 @@ END_DOC
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0 B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0
C_vector_DIIS(dim_DIIS+1) = -1.d0 C_vector_DIIS(dim_DIIS+1) = -1.d0
! Solve the linear system C = B.X deallocate(scratch)
! Estimate condition number of B
double precision :: anorm
integer :: info integer :: info
integer,allocatable :: ipiv(:) integer,allocatable :: ipiv(:)
allocate( &
ipiv(dim_DIIS+1) &
)
double precision, allocatable :: AF(:,:) double precision, allocatable :: AF(:,:)
allocate (AF(dim_DIIS+1,dim_DIIS+1)) double precision, external :: dlange
double precision :: rcond, ferr, berr
integer :: iwork(dim_DIIS+1), lwork
call dsysvx('N','U',dim_DIIS+1,1, & lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5)
B_matrix_DIIS,size(B_matrix_DIIS,1), & allocate(AF(dim_DIIS+1,dim_DIIS+1))
AF, size(AF,1), & allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) )
ipiv, &
C_vector_DIIS,size(C_vector_DIIS,1), &
X_vector_DIIS,size(X_vector_DIIS,1), &
rcond, &
ferr, &
berr, &
scratch,-1, &
iwork, &
info &
)
lwork = int(scratch(1,1))
deallocate(scratch)
allocate(scratch(lwork,1)) allocate(scratch(lwork,1))
anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, &
size(B_matrix_DIIS,1), scratch)
AF(:,:) = B_matrix_DIIS(:,:)
call dgetrf(dim_DIIS+1,dim_DIIS+1,AF,size(AF,1),ipiv,info)
if (info /= 0) then
dim_DIIS = 0
return
endif
call dgecon( '1', dim_DIIS+1, AF, &
size(AF,1), anorm, rcond, scratch, iwork, info )
if (info /= 0) then
dim_DIIS = 0
return
endif
if (rcond < 1.d-10) then
dim_DIIS = 0
return
endif
! Solve the linear system C = B.X
call dsysvx('N','U',dim_DIIS+1,1, & call dsysvx('N','U',dim_DIIS+1,1, &
B_matrix_DIIS,size(B_matrix_DIIS,1), & B_matrix_DIIS,size(B_matrix_DIIS,1), &
AF, size(AF,1), & AF, size(AF,1), &
ipiv, & ipiv, &
C_vector_DIIS,size(C_vector_DIIS,1), & C_vector_DIIS,size(C_vector_DIIS,1), &
X_vector_DIIS,size(X_vector_DIIS,1), & X_vector_DIIS,size(X_vector_DIIS,1), &
rcond, & rcond, &
ferr, & ferr, &
berr, & berr, &
scratch,size(scratch), & scratch,size(scratch), &
iwork, & iwork, &
info & info &
) )
deallocate(scratch,AF,iwork)
if(info < 0) then if(info < 0) then
stop 'bug in DIIS' stop 'bug in DIIS'
endif endif
if (rcond > 1.d-12) then
! Compute extrapolated Fock matrix ! Compute extrapolated Fock matrix
!$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200)
do j=1,ao_num do j=1,ao_num
do i=1,ao_num do i=1,ao_num
Fock_matrix_AO_(i,j) = 0.d0 Fock_matrix_AO_(i,j) = 0.d0
enddo enddo
do k=1,dim_DIIS do k=1,dim_DIIS
if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle
do i=1,ao_num do i=1,ao_num
Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + & Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + &
X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1) X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)
enddo
enddo
enddo enddo
!$OMP END PARALLEL DO enddo
enddo
else !$OMP END PARALLEL DO
dim_DIIS = 0
endif
end end