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

Fix pseudos NaN

Browse Source
This commit is contained in:
Anthony Scemama 2021-07-06 15:02:47 +02:00
parent 0eff94633d
commit ea92daeb45
3 changed files with 48 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
scripts/module/module_handler.py
View File

@ -31,7 +31,7 @@ try:
from docopt import docopt from docopt import docopt
from qp_path import QP_SRC, QP_ROOT, QP_PLUGINS, QP_EZFIO from qp_path import QP_SRC, QP_ROOT, QP_PLUGINS, QP_EZFIO
except ImportError: except ImportError:
print("source .quantum_package.rc") print("source quantum_package.rc")
raise raise

59
src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
View File

@ -37,26 +37,58 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
integer :: num_A,num_B integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3) integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m integer :: i,j,k,l,m
double precision :: Vloc, Vpseudo double precision :: Vloc, Vpseudo
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 double precision :: wall_1, wall_2, wall_0
integer :: thread_num integer :: thread_num
integer :: omp_get_thread_num integer :: omp_get_thread_num
double precision :: c
double precision :: Z
PROVIDE ao_coef_normalized_ordered_transp
PROVIDE pseudo_v_k_transp pseudo_n_k_transp pseudo_klocmax pseudo_dz_k_transp
ao_pseudo_integrals_local = 0.d0 ao_pseudo_integrals_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (local)' print*, 'Providing the nuclear electron pseudo integrals (local)'
call wall_time(wall_1) ! Dummy iteration for OpenMP
call cpu_time(cpu_1) j=1
i=1
l=1
m=1
num_A = ao_nucl(j)
power_A(1:3)= ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(i)
power_B(1:3)= ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
alpha = ao_expo_ordered_transp(l,j)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
c = c + Vloc(pseudo_klocmax, &
pseudo_v_k_transp (1,k), &
pseudo_n_k_transp (1,k), &
pseudo_dz_k_transp(1,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
ao_pseudo_integrals_local = 0.d0
call wall_time(wall_1)
thread_num = 0 thread_num = 0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
!$OMP num_A,num_B,Z,c,n_pt_in, & !$OMP num_A,num_B,Z,c, &
!$OMP wall_0,wall_2,thread_num) & !$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP ao_pseudo_integrals_local,nucl_num,nucl_charge, & !$OMP ao_pseudo_integrals_local,nucl_num,nucl_charge, &
@ -66,7 +98,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
!$ thread_num = omp_get_thread_num() !$ thread_num = omp_get_thread_num()
wall_0 = wall_1 wall_0 = wall_1
!$OMP DO SCHEDULE (guided) !$OMP DO
do j = 1, ao_num do j = 1, ao_num
@ -85,7 +117,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
do m=1,ao_prim_num(i) do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i) beta = ao_expo_ordered_transp(m,i)
double precision :: c
c = 0.d0 c = 0.d0
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
@ -93,7 +124,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
cycle cycle
endif endif
do k = 1, nucl_num do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k) Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)
@ -137,25 +167,28 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
integer :: num_A,num_B integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3) integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m integer :: i,j,k,l,m
double precision :: Vloc, Vpseudo double precision :: Vloc, Vpseudo
integer :: omp_get_thread_num integer :: omp_get_thread_num
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 double precision :: wall_1, wall_2, wall_0
integer :: thread_num integer :: thread_num
double precision :: c
double precision :: Z
PROVIDE ao_coef_normalized_ordered_transp
PROVIDE pseudo_lmax pseudo_kmax pseudo_v_kl_transp pseudo_n_kl_transp pseudo_dz_kl_transp
ao_pseudo_integrals_non_local = 0.d0 ao_pseudo_integrals_non_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (non-local)' print*, 'Providing the nuclear electron pseudo integrals (non-local)'
call wall_time(wall_1) call wall_time(wall_1)
call cpu_time(cpu_1)
thread_num = 0 thread_num = 0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
!$OMP num_A,num_B,Z,c,n_pt_in, & !$OMP num_A,num_B,Z,c, &
!$OMP wall_0,wall_2,thread_num) & !$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP ao_pseudo_integrals_non_local,nucl_num,nucl_charge,& !$OMP ao_pseudo_integrals_non_local,nucl_num,nucl_charge,&
@ -184,7 +217,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
do m=1,ao_prim_num(i) do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i) beta = ao_expo_ordered_transp(m,i)
double precision :: c
c = 0.d0 c = 0.d0
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
@ -193,7 +225,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
endif endif
do k = 1, nucl_num do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k) Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)

4
src/ao_one_e_ints/pseudopot.f90
View File

@ -666,7 +666,7 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2)
arg=g_a*ac**2+g_b*bc**2 arg=g_a*ac**2+g_b*bc**2
if(arg.gt.-dlog(10.d-20))then if(arg.gt.-dlog(1.d-20))then
Vloc=0.d0 Vloc=0.d0
return return
endif endif
@ -1839,7 +1839,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
m_1 = m+m+1 m_1 = m+m+1
nlm = n+m+l nlm = n+m+l
pi=dacos(-1.d0) pi=dacos(-1.d0)
a_over_b_square = (a/b)**2 a_over_b_square = (a*a)/(b*b)
! First term of the sequence ! First term of the sequence