From a90b446bebfb3410fef76fe2715819197fa2d16d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 May 2020 18:11:27 +0200 Subject: [PATCH] Fixed floating-point exceptions --- bin/qp_reset | 2 +- config/gfortran_debug.cfg | 2 +- config/travis.cfg | 2 +- src/ao_one_e_ints/pseudopot.f90 | 6 +- src/hartree_fock/10.hf.bats | 13 +-- src/nuclei/nuclei.irp.f | 24 +++--- src/scf_utils/roothaan_hall_scf.irp.f | 120 ++++++++++++++------------ 7 files changed, 93 insertions(+), 76 deletions(-) diff --git a/bin/qp_reset b/bin/qp_reset index 70cc7c9f..74dd1f78 100755 --- a/bin/qp_reset +++ b/bin/qp_reset @@ -112,7 +112,7 @@ qp_edit --check ${ezfio} if [[ $mos -eq 1 ]] ; then 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_pseudo None qp set mo_one_e_ints io_mo_one_e_integrals None diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 0acb5fa5..342acae9 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -51,7 +51,7 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [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 ################# diff --git a/config/travis.cfg b/config/travis.cfg index 93e63f7b..ca789238 100644 --- a/config/travis.cfg +++ b/config/travis.cfg @@ -53,7 +53,7 @@ FCFLAGS : -Ofast -fimplicit-none # -g : Extra debugging information # [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 diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 index 563705dd..31c3c549 100644 --- a/src/ao_one_e_ints/pseudopot.f90 +++ b/src/ao_one_e_ints/pseudopot.f90 @@ -1907,13 +1907,17 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) sum=s_q_0 if (q>300) then - stop 'pseudopot.f90 : q > 300' + stop 'pseudopot.f90 : q > 200' endif qk = dble(q) two_qkmp1 = 2.d0*(qk+mk)+1.d0 do k=0,q-1 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 two_qkmp1 = two_qkmp1-2.d0 qk = qk-1.d0 diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index 0301748b..e3e000a9 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -5,12 +5,13 @@ source $QP_ROOT/quantum_package.rc function run() { - thresh=1.e-7 + thresh=1.e-6 test_exe scf || skip qp set_file $1 qp edit --check qp reset --mos - qp run scf + qp set scf_utils n_it_scf_max 50 + qp run scf # qp set_frozen_core energy="$(ezfio get hartree_fock energy)" eq $energy $2 $thresh @@ -39,7 +40,7 @@ function run() { } @test "SO" { # 0.539000 5.70403s - run so.ezfio -25.7175263371941 + run so.ezfio -25.7175270084056 } @test "HCO" { # 0.636700 1.55279s @@ -107,13 +108,13 @@ function run() { } @test "C2H2" { # 19.599000 37.7923s - run c2h2.ezfio -12.12144019495306 + run c2h2.ezfio -12.12144044853196 } @test "SiH3" { # 20.316100 54.0861s [[ -n $TRAVIS ]] && skip - run sih3.ezfio -5.455398769158780 + run sih3.ezfio -5.455400439077580 } @test "OH" { # 32.042200 1.36478m @@ -130,6 +131,6 @@ function run() { @test "SO2" { # 71.894900 3.22567m [[ -n $TRAVIS ]] && skip - run so2.ezfio -41.55800190733211 + run so2.ezfio -41.55800401346361 } diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index 523bd80f..c1b5f52f 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -53,18 +53,20 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] '================','============','============','============','============' write(6,'(A)') '' - double precision :: dist_min, x, y, z - dist_min = huge(1.d0) - do i=1,nucl_num - do j=i+1,nucl_num - x = nucl_coord(i,1)-nucl_coord(j,1) - y = nucl_coord(i,2)-nucl_coord(j,2) - z = nucl_coord(i,3)-nucl_coord(j,3) - dist_min = min(x*x + y*y + z*z, dist_min) + if (nucl_num > 1) then + double precision :: dist_min, x, y, z + dist_min = huge(1.d0) + do i=1,nucl_num + do j=i+1,nucl_num + x = nucl_coord(i,1)-nucl_coord(j,1) + y = nucl_coord(i,2)-nucl_coord(j,2) + z = nucl_coord(i,3)-nucl_coord(j,3) + dist_min = min(x*x + y*y + z*z, dist_min) + enddo enddo - enddo - write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', & - dsqrt(dist_min)*a0,' Angstrom' + write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', & + dsqrt(dist_min)*a0,' Angstrom' + endif endif diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index ea472cdf..8c76501e 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -196,6 +196,13 @@ END_DOC double precision,allocatable :: scratch(:,:) 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( & 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 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,allocatable :: ipiv(:) - - allocate( & - ipiv(dim_DIIS+1) & - ) - double precision, allocatable :: AF(:,:) - allocate (AF(dim_DIIS+1,dim_DIIS+1)) - double precision :: rcond, ferr, berr - integer :: iwork(dim_DIIS+1), lwork + double precision, external :: dlange - call dsysvx('N','U',dim_DIIS+1,1, & - B_matrix_DIIS,size(B_matrix_DIIS,1), & - AF, size(AF,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) + lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5) + allocate(AF(dim_DIIS+1,dim_DIIS+1)) + allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+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, & - B_matrix_DIIS,size(B_matrix_DIIS,1), & - AF, size(AF,1), & - ipiv, & - C_vector_DIIS,size(C_vector_DIIS,1), & - X_vector_DIIS,size(X_vector_DIIS,1), & - rcond, & - ferr, & - berr, & - scratch,size(scratch), & - iwork, & - info & - ) + B_matrix_DIIS,size(B_matrix_DIIS,1), & + AF, size(AF,1), & + ipiv, & + C_vector_DIIS,size(C_vector_DIIS,1), & + X_vector_DIIS,size(X_vector_DIIS,1), & + rcond, & + ferr, & + berr, & + scratch,size(scratch), & + iwork, & + info & + ) + + deallocate(scratch,AF,iwork) if(info < 0) then stop 'bug in DIIS' endif - if (rcond > 1.d-12) then - ! Compute extrapolated Fock matrix - !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) - do j=1,ao_num - do i=1,ao_num - Fock_matrix_AO_(i,j) = 0.d0 - enddo - do k=1,dim_DIIS - if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle - do i=1,ao_num - Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + & - X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1) - enddo - enddo + !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) + do j=1,ao_num + do i=1,ao_num + Fock_matrix_AO_(i,j) = 0.d0 + enddo + do k=1,dim_DIIS + if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle + do i=1,ao_num + Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + & + X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1) enddo - !$OMP END PARALLEL DO - - else - dim_DIIS = 0 - endif + enddo + enddo + !$OMP END PARALLEL DO end