From bf2e02dc9dbc4aa482336f1b25f577fc4de2b925 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sun, 16 Apr 2017 22:49:38 +0200 Subject: [PATCH] beginning to clean erf integrals --- plugins/Integrals_erf/EZFIO.cfg | 16 ++----- .../map_integrals_long_range.irp.f | 19 +++++++- plugins/Integrals_erf/mo_bi_integrals.irp.f | 6 +-- .../Integrals_erf/mo_bi_integrals_erf.irp.f | 30 ++++++------ plugins/Integrals_erf/providers_ao_erf.irp.f | 32 ++++++------- .../Integrals_erf/providers_ao_standard.irp.f | 2 +- plugins/Integrals_erf/read_write.irp.f | 48 +++++++++++++++++++ src/Determinants/EZFIO.cfg | 6 +++ src/Determinants/density_matrix.irp.f | 15 ++++++ src/Integrals_Bielec/EZFIO.cfg | 15 ------ src/Integrals_Bielec/ao_bi_integrals.irp.f | 37 ++++---------- src/Integrals_Monoelec/mo_mono_ints.irp.f | 25 +++++++--- 12 files changed, 151 insertions(+), 100 deletions(-) diff --git a/plugins/Integrals_erf/EZFIO.cfg b/plugins/Integrals_erf/EZFIO.cfg index fb08245f..ca798f68 100644 --- a/plugins/Integrals_erf/EZFIO.cfg +++ b/plugins/Integrals_erf/EZFIO.cfg @@ -32,16 +32,16 @@ doc: Read/Write MO integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[disk_access_ao_integrals_standard] +[disk_access_ao_integrals_erf] type: Disk_access -doc: Read/Write AO integrals_standard from/to disk [ Write | Read | None ] +doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[disk_access_mo_integrals_standard] +[disk_access_mo_integrals_erf] type: Disk_access -doc: Read/Write MO integrals_standard from/to disk [ Write | Read | None ] +doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None @@ -74,11 +74,3 @@ interface: ezfio,provider,ocaml default: 0.5 ezfio_name: mu_erf - -[long_range] -type: logical -doc: if true, compute all the integrals using the long range interaction -interface: ezfio,provider,ocaml -default: False -ezfio_name: long_range - diff --git a/plugins/Integrals_erf/map_integrals_long_range.irp.f b/plugins/Integrals_erf/map_integrals_long_range.irp.f index 2ab2f14f..8f24692a 100644 --- a/plugins/Integrals_erf/map_integrals_long_range.irp.f +++ b/plugins/Integrals_erf/map_integrals_long_range.irp.f @@ -28,6 +28,7 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ] + use map_module implicit none BEGIN_DOC ! Cache of AO integrals for fast access @@ -42,7 +43,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ] do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max !DIR$ FORCEINLINE - call bielec_integrals_erf_index(i,j,k,l,idx) + call bielec_integrals_index(i,j,k,l,idx) !DIR$ FORCEINLINE call map_get(ao_integrals_erf_map,idx,integral) ii = l-ao_integrals_erf_cache_min @@ -343,7 +344,21 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ] print*, 'MO map initialized' END_PROVIDER -subroutine insert_into_mo_integrals_map(n_integrals, & +subroutine insert_into_ao_integrals_erf_map(n_integrals,buffer_i, buffer_values) + use map_module + implicit none + BEGIN_DOC + ! Create new entry into AO map + END_DOC + + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + + call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals) +end + +subroutine insert_into_mo_integrals_erf_map(n_integrals, & buffer_i, buffer_values, thr) use map_module implicit none diff --git a/plugins/Integrals_erf/mo_bi_integrals.irp.f b/plugins/Integrals_erf/mo_bi_integrals.irp.f index f4a20cfc..9cca2e12 100644 --- a/plugins/Integrals_erf/mo_bi_integrals.irp.f +++ b/plugins/Integrals_erf/mo_bi_integrals.irp.f @@ -122,7 +122,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] if (write_mo_integrals) then call ezfio_set_work_empty(.False.) call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) - call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + call ezfio_set_integrals_erf_disk_access_mo_integrals("Read") endif END_PROVIDER @@ -1065,8 +1065,6 @@ end ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij END_DOC implicit none - long_range = .false. - touch long_range ! call routine_mo_bielec_integral_jj_from_ao ! END_PROVIDER @@ -1199,8 +1197,6 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] &BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] implicit none - long_range = .false. - touch long_range BEGIN_DOC ! mo_bielec_integral_vv_from_ao(i,j) = J_ij ! mo_bielec_integral_vv_exchange_from_ao(i,j) = J_ij diff --git a/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f b/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f index c7193819..ecb55757 100644 --- a/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f +++ b/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f @@ -30,14 +30,14 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ] END_DOC mo_bielec_integrals_erf_in_map = .True. -! if (read_mo_integrals) then -! print*,'Reading the MO integrals' -! call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) -! print*, 'MO integrals provided' -! return -! else - PROVIDE ao_bielec_integrals_in_map -! endif + if (read_mo_integrals_erf) then + print*,'Reading the MO integrals_erf' + call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map) + print*, 'MO integrals_erf provided' + return + else + PROVIDE ao_bielec_integrals_erf_in_map + endif !if(no_vvvv_integrals)then ! integer :: i,j,k,l @@ -117,13 +117,13 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ] ! endif ! !else - call add_integrals_to_map(full_ijkl_bitmask_4) - !endif - !if (write_mo_integrals) then - ! call ezfio_set_work_empty(.False.) - ! call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) - ! call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + call add_integrals_erf_to_map(full_ijkl_bitmask_4) !endif + if (write_mo_integrals_erf) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map) + call ezfio_set_integrals_erf_disk_access_mo_integrals_erf("Read") + endif END_PROVIDER @@ -160,7 +160,7 @@ subroutine add_integrals_erf_to_map(mask_ijkl) integer :: i2,i3,i4 double precision,parameter :: thr_coef = 1.d-10 - PROVIDE ao_bielec_integrals_in_map mo_coef + PROVIDE ao_bielec_integrals_erf_in_map mo_coef !Get list of MOs for i,j,k and l !------------------------------- diff --git a/plugins/Integrals_erf/providers_ao_erf.irp.f b/plugins/Integrals_erf/providers_ao_erf.irp.f index cff3e49c..1507d1be 100644 --- a/plugins/Integrals_erf/providers_ao_erf.irp.f +++ b/plugins/Integrals_erf/providers_ao_erf.irp.f @@ -25,22 +25,22 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ] integral = ao_bielec_integral_erf(1,1,1,1) real :: map_mb -! PROVIDE read_ao_integrals disk_access_ao_integrals -! if (read_ao_integrals) then -! print*,'Reading the AO integrals' -! call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) -! print*, 'AO integrals provided' -! ao_bielec_integrals_in_map = .True. -! return -! endif + PROVIDE read_ao_integrals_erf disk_access_ao_integrals_erf + if (read_ao_integrals_erf) then + print*,'Reading the AO integrals_erf' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map) + print*, 'AO integrals_erf provided' + ao_bielec_integrals_erf_in_map = .True. + return + endif - print*, 'Providing the AO integrals' + print*, 'Providing the AO integrals_erf' call wall_time(wall_0) call wall_time(wall_1) call cpu_time(cpu_1) integer(ZMQ_PTR) :: zmq_to_qp_run_socket - call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') + call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals_erf') character(len=:), allocatable :: task allocate(character(len=ao_num*12) :: task) @@ -63,7 +63,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ] endif !$OMP END PARALLEL - call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals') + call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals_erf') print*, 'Sorting the map' @@ -81,11 +81,11 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ] ao_bielec_integrals_erf_in_map = .True. -! if (write_ao_integrals) then -! call ezfio_set_work_empty(.False.) -! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map) -! call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read") -! endif + if (write_ao_integrals_erf) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map) + call ezfio_set_integrals_erf_disk_access_ao_integrals_erf("Read") + endif END_PROVIDER diff --git a/plugins/Integrals_erf/providers_ao_standard.irp.f b/plugins/Integrals_erf/providers_ao_standard.irp.f index c0b043da..973c0e58 100644 --- a/plugins/Integrals_erf/providers_ao_standard.irp.f +++ b/plugins/Integrals_erf/providers_ao_standard.irp.f @@ -84,7 +84,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] if (write_ao_integrals) then call ezfio_set_work_empty(.False.) call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read") + call ezfio_set_integrals_erf_disk_access_ao_integrals("Read") endif END_PROVIDER diff --git a/plugins/Integrals_erf/read_write.irp.f b/plugins/Integrals_erf/read_write.irp.f index 5b2b7f3e..e9fc0f91 100644 --- a/plugins/Integrals_erf/read_write.irp.f +++ b/plugins/Integrals_erf/read_write.irp.f @@ -45,3 +45,51 @@ implicit none endif END_PROVIDER + +BEGIN_PROVIDER [ logical, read_ao_integrals_erf ] +&BEGIN_PROVIDER [ logical, read_mo_integrals_erf ] +&BEGIN_PROVIDER [ logical, write_ao_integrals_erf ] +&BEGIN_PROVIDER [ logical, write_mo_integrals_erf ] + + BEGIN_DOC +! One level of abstraction for disk_access_ao_integrals_erf and disk_access_mo_integrals_erf + END_DOC +implicit none + + if (disk_access_ao_integrals_erf.EQ.'Read') then + read_ao_integrals_erf = .True. + write_ao_integrals_erf = .False. + + else if (disk_access_ao_integrals_erf.EQ.'Write') then + read_ao_integrals_erf = .False. + write_ao_integrals_erf = .True. + + else if (disk_access_ao_integrals_erf.EQ.'None') then + read_ao_integrals_erf = .False. + write_ao_integrals_erf = .False. + + else + print *, 'bielec_integrals_erf/disk_access_ao_integrals_erf has a wrong type' + stop 1 + + endif + + if (disk_access_mo_integrals_erf.EQ.'Read') then + read_mo_integrals_erf = .True. + write_mo_integrals_erf = .False. + + else if (disk_access_mo_integrals_erf.EQ.'Write') then + read_mo_integrals_erf = .False. + write_mo_integrals_erf = .True. + + else if (disk_access_mo_integrals_erf.EQ.'None') then + read_mo_integrals_erf = .False. + write_mo_integrals_erf = .False. + + else + print *, 'bielec_integrals_erf/disk_access_mo_integrals_erf has a wrong type' + stop 1 + + endif + +END_PROVIDER diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index a68a61a5..a9ecd806 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -119,3 +119,9 @@ doc: Maximum number of determinants for which the full H matrix is stored. Be ca interface: ezfio,provider,ocaml default: 90000 +[density_matrix_mo_disk] +type: double precision +doc: coefficient of the ith ao on the jth mo +interface: ezfio +size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num) + diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 7274760e..541cfcb4 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -15,6 +15,21 @@ enddo END_PROVIDER + +subroutine save_density_matrix_mo + implicit none + double precision, allocatable :: dm(:,:) + allocate(dm(mo_tot_num,mo_tot_num)) + integer :: i,j + do i = 1, mo_tot_num + do j = 1, mo_tot_num + dm(i,j) = one_body_dm_mo_alpha_average(i,j) + enddo + enddo + call ezfio_set_determinants_density_matrix_mo_disk(dm) + +end + BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] implicit none integer :: i,j,ispin,istate diff --git a/src/Integrals_Bielec/EZFIO.cfg b/src/Integrals_Bielec/EZFIO.cfg index db965d69..0576b811 100644 --- a/src/Integrals_Bielec/EZFIO.cfg +++ b/src/Integrals_Bielec/EZFIO.cfg @@ -52,18 +52,3 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_mo -[mu_erf] -type: double precision -doc: cutting of the interaction in the range separated model -interface: ezfio,provider,ocaml -default: 0.5 -ezfio_name: mu_erf - - -[long_range] -type: logical -doc: if true, compute all the integrals using the long range interaction -interface: ezfio,provider,ocaml -default: False -ezfio_name: long_range - diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 54719612..196bfce4 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -158,7 +158,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) q_inv = 1.d0/qq schwartz_kl(s,r) = general_primitive_integral(dim1, & Q_new,Q_center,fact_q,qq,q_inv,iorder_q, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) & * coef2 schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) enddo @@ -471,15 +471,9 @@ double precision function general_primitive_integral(dim, & ! Gaussian Product ! ---------------- - double precision :: p_plus_q - if(long_range)then - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - else - p_plus_q = (p+q) - endif - pq = p_inv*0.5d0*q_inv - pq_inv = 0.5d0/p_plus_q + pq = p_inv*0.5d0*q_inv + pq_inv = 0.5d0/(p+q) p10_1 = q*pq ! 1/(2p) p01_1 = p*pq ! 1/(2q) pq_inv_2 = pq_inv+pq_inv @@ -554,7 +548,7 @@ double precision function general_primitive_integral(dim, & return endif - rho = p*q *pq_inv_2 ! le rho qui va bien + rho = p*q *pq_inv_2 dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + & (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + & (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3)) @@ -580,8 +574,7 @@ double precision function general_primitive_integral(dim, & double precision :: rint_sum accu = accu + rint_sum(n_pt_out,const,d1) - ! change p+q in dsqrt - general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q) + general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q) end @@ -627,16 +620,10 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y, ASSERT (gama >= 0.d0) p = alpha + beta q = delta + gama - double precision :: p_plus_q - if(long_range)then - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - else - p_plus_q = (p+q) - endif ASSERT (p+q >= 0.d0) n_pt = ishft( nx+ny+nz,1 ) - coeff = pi_5_2 / (p * q * dsqrt(p_plus_q)) + coeff = pi_5_2 / (p * q * dsqrt(p+q)) if (n_pt == 0) then ERI = coeff return @@ -663,21 +650,15 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q integer :: i, n_iter, n_pt, j double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2 integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz - double precision :: p_plus_q j = ishft(n_pt,-1) ASSERT (n_pt > 1) - if(long_range)then - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - else - p_plus_q = (p+q) - endif - pq_inv = 0.5d0/(p_plus_q) + pq_inv = 0.5d0/(p+q) pq_inv_2 = pq_inv + pq_inv p10_1 = 0.5d0/p p01_1 = 0.5d0/q - p10_2 = 0.5d0 * q /(p * p_plus_q) - p01_2 = 0.5d0 * p /(q * p_plus_q) + p10_2 = 0.5d0 * q /(p * q + p * p) + p01_2 = 0.5d0 * p /(q * q + q * p) double precision :: B00(n_pt_max_integrals) double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals) double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) diff --git a/src/Integrals_Monoelec/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f index 50ab7ffa..20aae1fd 100644 --- a/src/Integrals_Monoelec/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -6,10 +6,23 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to ! sum of the kinetic and nuclear electronic potential END_DOC print*,'Providing the mono electronic integrals' - do j = 1, mo_tot_num - do i = 1, mo_tot_num - mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + & - mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j) - enddo - enddo + if (read_mo_one_integrals) then + call read_one_e_integrals('mo_one_integral', mo_mono_elec_integral, & + size(mo_mono_elec_integral,1), size(mo_mono_elec_integral,2)) + print *, 'MO N-e integrals read from disk' + else + do j = 1, mo_tot_num + do i = 1, mo_tot_num + mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + & + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j) + enddo + enddo + endif + +! if (write_mo_one_integrals) then +! call write_one_e_integrals('mo_one_integral', mo_mono_elec_integral, & +! size(mo_mono_elec_integral,1), size(mo_mono_elec_integral,2)) +! print *, 'MO N-e integrals written to disk' +! endif + END_PROVIDER