From bb23d6a5b5160387ef0695d1a07e7f4ef86f71b6 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 12 Jun 2023 13:36:01 +0200 Subject: [PATCH 1/4] Fixed the pt_charges bug: + added the pt_charges integrals to the usual v_ne + added only the nuclei_pt_charge interaction to the usual nuclear_repulsion (and not the pt_charge_pt_charge interaction) --- src/ao_one_e_ints/pot_ao_ints.irp.f | 3 +++ src/hartree_fock/10.hf.bats | 5 +++-- src/nuclei/nuclei.irp.f | 7 ++++++- src/nuclei/point_charges.irp.f | 3 +++ 4 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 446bf730..4f9ae76d 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -104,6 +104,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] IF(do_pseudo) THEN ao_integrals_n_e += ao_pseudo_integrals ENDIF + IF(point_charges) THEN + ao_integrals_n_e += ao_integrals_pt_chrg + ENDIF endif diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index 3647b775..6e7d0233 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -43,12 +43,11 @@ python write_pt_charges.py ${EZFIO} qp set nuclei point_charges True qp run scf | tee ${EZFIO}.pt_charges.out energy="$(ezfio get hartree_fock energy)" -good=-92.76613324421798 +good=-92.79920682236470 eq $energy $good $thresh rm -rf $EZFIO } - @test "H2_1" { # 1s run h2_1.ezfio -1.005924963288527 } @@ -85,6 +84,8 @@ rm -rf $EZFIO run hcn.ezfio -92.88717500035233 } + + @test "B-B" { # 3s run b2_stretched.ezfio -48.9950585434279 } diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index fabdc42e..bb8cc782 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -206,7 +206,12 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] enddo nuclear_repulsion *= 0.5d0 if(point_charges)then - nuclear_repulsion += pt_chrg_nuclei_interaction + pt_chrg_interaction + print*,'bear nuclear repulsion = ',nuclear_repulsion + print*,'adding the interaction between the nuclein and the point charges' + print*,'to the usual nuclear repulsion ' + nuclear_repulsion += pt_chrg_nuclei_interaction + print*,'new nuclear repulsion = ',nuclear_repulsion + print*,'WARNING: we do not add the interaction between the point charges themselves' endif end if diff --git a/src/nuclei/point_charges.irp.f b/src/nuclei/point_charges.irp.f index b955537f..66905c8c 100644 --- a/src/nuclei/point_charges.irp.f +++ b/src/nuclei/point_charges.irp.f @@ -205,5 +205,8 @@ BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_interaction] enddo print*,'Interaction between point charges and nuclei' print*,'pt_chrg_nuclei_interaction = ',pt_chrg_nuclei_interaction + if(point_charges)then + provide pt_chrg_interaction + endif END_PROVIDER From 4d9e28438c199c7f8956913b7380c5ba6ec07932 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 12 Jun 2023 14:05:36 +0200 Subject: [PATCH 2/4] Improved I/O in CCSD --- src/ao_two_e_ints/cholesky.irp.f | 13 ++---- src/ccsd/EZFIO.cfg | 11 +++++ src/ccsd/ccsd_space_orb_sub.irp.f | 11 +++-- src/ccsd/ccsd_spin_orb_sub.irp.f | 14 +++--- src/ccsd/save_energy.irp.f | 13 ++++++ src/mo_two_e_ints/cholesky.irp.f | 2 + src/utils/linear_algebra.irp.f | 6 +-- src/utils_cc/EZFIO.cfg | 16 +++---- src/utils_cc/guess_t.irp.f | 75 +++++++++++++++---------------- 9 files changed, 88 insertions(+), 73 deletions(-) create mode 100644 src/ccsd/EZFIO.cfg create mode 100644 src/ccsd/save_energy.irp.f diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index bb81b141..77eb6ddc 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ] ! Number of Cholesky vectors in AO basis END_DOC - cholesky_ao_num_guess = ao_num*ao_num / 2 + cholesky_ao_num_guess = ao_num*ao_num END_PROVIDER BEGIN_PROVIDER [ integer, cholesky_ao_num ] @@ -44,19 +44,12 @@ END_PROVIDER do m=0,9 do l=1+m,ao_num,10 !$OMP DO SCHEDULE(dynamic) - do j=1,l + do j=1,ao_num do k=1,ao_num - do i=1,min(k,j) + do i=1,ao_num if (ao_two_e_integral_zero(i,j,k,l)) cycle integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map) ao_integrals(i,k,j,l) = integral - ao_integrals(k,i,j,l) = integral - ao_integrals(i,k,l,j) = integral - ao_integrals(k,i,l,j) = integral - ao_integrals(j,l,i,k) = integral - ao_integrals(j,l,k,i) = integral - ao_integrals(l,j,i,k) = integral - ao_integrals(l,j,k,i) = integral enddo enddo enddo diff --git a/src/ccsd/EZFIO.cfg b/src/ccsd/EZFIO.cfg new file mode 100644 index 00000000..328cd981 --- /dev/null +++ b/src/ccsd/EZFIO.cfg @@ -0,0 +1,11 @@ +[energy] +type: double precision +doc: CCSD energy +interface: ezfio + +[energy_t] +type: double precision +doc: CCSD(T) energy +interface: ezfio + + diff --git a/src/ccsd/ccsd_space_orb_sub.irp.f b/src/ccsd/ccsd_space_orb_sub.irp.f index 1467d9a4..40c57188 100644 --- a/src/ccsd/ccsd_space_orb_sub.irp.f +++ b/src/ccsd/ccsd_space_orb_sub.irp.f @@ -135,8 +135,11 @@ subroutine run_ccsd_space_orb write(*,'(A15,1pE10.2,A3)')' Conv = ', max_r print*,'' - call write_t1(nO,nV,t1) - call write_t2(nO,nV,t2) + if (write_amplitudes) then + call write_t1(nO,nV,t1) + call write_t2(nO,nV,t2) + call ezfio_set_utils_cc_io_amplitudes('Read') + endif ! Deallocation if (cc_update_method == 'diis') then @@ -147,6 +150,7 @@ subroutine run_ccsd_space_orb ! CCSD(T) double precision :: e_t + e_t = 0.d0 if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then @@ -182,8 +186,7 @@ subroutine run_ccsd_space_orb print*,'' endif - print*,'Reference determinant:' - call print_det(det,N_int) + call save_energy(uncorr_energy + energy, e_t) deallocate(t1,t2) diff --git a/src/ccsd/ccsd_spin_orb_sub.irp.f b/src/ccsd/ccsd_spin_orb_sub.irp.f index 23e2cef1..a267cc45 100644 --- a/src/ccsd/ccsd_spin_orb_sub.irp.f +++ b/src/ccsd/ccsd_spin_orb_sub.irp.f @@ -269,8 +269,11 @@ subroutine run_ccsd_spin_orb write(*,'(A15,1pE10.2,A3)')' Conv = ', max_r print*,'' - call write_t1(nO,nV,t1) - call write_t2(nO,nV,t2) + if (write_amplitudes) then + call write_t1(nO,nV,t1) + call write_t2(nO,nV,t2) + call ezfio_set_utils_cc_io_amplitudes('Read') + endif ! Deallocate if (cc_update_method == 'diis') then @@ -284,8 +287,9 @@ subroutine run_ccsd_spin_orb deallocate(v_ovoo,v_oovo) deallocate(v_ovvo,v_ovov,v_oovv) + double precision :: t_corr + t_corr = 0.d0 if (cc_par_t .and. elec_alpha_num +elec_beta_num > 2) then - double precision :: t_corr print*,'CCSD(T) calculation...' call wall_time(ta) !allocate(v_vvvo(nV,nV,nV,nO)) @@ -307,8 +311,8 @@ subroutine run_ccsd_spin_orb write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + t_corr, ' Ha' print*,'' endif - print*,'Reference determinant:' - call print_det(det,N_int) + + call save_energy(uncorr_energy + energy, t_corr) deallocate(f_oo,f_ov,f_vv,f_o,f_v) deallocate(v_ooov,v_vvoo,t1,t2) diff --git a/src/ccsd/save_energy.irp.f b/src/ccsd/save_energy.irp.f new file mode 100644 index 00000000..30d93ec3 --- /dev/null +++ b/src/ccsd/save_energy.irp.f @@ -0,0 +1,13 @@ +subroutine save_energy(E,ET) + implicit none + BEGIN_DOC +! Saves the energy in |EZFIO|. + END_DOC + double precision, intent(in) :: E, ET + call ezfio_set_ccsd_energy(E) + if (ET /= 0.d0) then + call ezfio_set_ccsd_energy_t(E+ET) + endif +end + + diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f index 8b1e6e1c..32c0dccd 100644 --- a/src/mo_two_e_ints/cholesky.irp.f +++ b/src/mo_two_e_ints/cholesky.irp.f @@ -27,6 +27,8 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_ao_num, mo_num, double precision, allocatable :: buffer(:,:) print *, 'AO->MO Transformation of Cholesky vectors .' + + call set_multiple_levels_omp(.False.) !$OMP PARALLEL PRIVATE(i,j,k,buffer) allocate(buffer(mo_num,mo_num)) !$OMP DO SCHEDULE(static) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 69873bc0..76a539a6 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1831,7 +1831,7 @@ double precision, intent(in) :: tol integer, dimension(:), allocatable :: piv double precision, dimension(:), allocatable :: work -character, parameter :: uplo = "U" +character, parameter :: uplo = 'L' integer :: LDA integer :: info integer :: k, l, rank0 @@ -1848,14 +1848,14 @@ if (rank > rank0) then end if do k = 1, ndim - A(k+1:ndim, k) = 0.00D+0 + A(k,k+1:ndim) = 0.00D+0 end do ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer ! to do the swapping in-place U(:,:) = 0.00D+0 do k = 1, ndim l = piv(k) - U(l, 1:rank) = A(1:rank, k) + U(l, 1:rank) = A(k,1:rank) end do end subroutine pivoted_cholesky diff --git a/src/utils_cc/EZFIO.cfg b/src/utils_cc/EZFIO.cfg index 71ee87e3..fb6d9034 100644 --- a/src/utils_cc/EZFIO.cfg +++ b/src/utils_cc/EZFIO.cfg @@ -46,17 +46,11 @@ doc: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation t interface: ezfio,ocaml,provider default: MP -[cc_write_t1] -type: logical -doc: If true, it will write on disk the T1 amplitudes at the end of the calculation. -interface: ezfio,ocaml,provider -default: False - -[cc_write_t2] -type: logical -doc: If true, it will write on disk the T2 amplitudes at the end of the calculation. -interface: ezfio,ocaml,provider -default: False +[io_amplitudes] +type: Disk_access +doc: Read/Write |CCSD| amplitudes from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None [cc_par_t] type: logical diff --git a/src/utils_cc/guess_t.irp.f b/src/utils_cc/guess_t.irp.f index 42acdf78..bb26e133 100644 --- a/src/utils_cc/guess_t.irp.f +++ b/src/utils_cc/guess_t.irp.f @@ -91,16 +91,17 @@ subroutine write_t1(nO,nV,t1) double precision, intent(in) :: t1(nO, nV) ! internal - integer :: i,a + integer :: i,a, iunit + integer, external :: getunitandopen - if (cc_write_t1) then - open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1') + if (write_amplitudes) then + iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','w') do a = 1, nV do i = 1, nO - write(11,'(F20.12)') t1(i,a) + write(iunit,'(F20.12)') t1(i,a) enddo enddo - close(11) + close(iunit) endif end @@ -120,20 +121,21 @@ subroutine write_t2(nO,nV,t2) double precision, intent(in) :: t2(nO, nO, nV, nV) ! internal - integer :: i,j,a,b + integer :: i,j,a,b, iunit + integer, external :: getunitandopen - if (cc_write_t2) then - open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2') + if (write_amplitudes) then + iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','w') do b = 1, nV do a = 1, nV do j = 1, nO do i = 1, nO - write(11,'(F20.12)') t2(i,j,a,b) + write(iunit,'(F20.12)') t2(i,j,a,b) enddo enddo enddo enddo - close(11) + close(iunit) endif end @@ -153,23 +155,19 @@ subroutine read_t1(nO,nV,t1) double precision, intent(out) :: t1(nO, nV) ! internal - integer :: i,a + integer :: i,a, iunit logical :: ok + integer, external :: getunitandopen - inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok) - if (.not. ok) then - print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1' - print*, 'Do a first calculation with cc_write_t1 = True' - print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read' - call abort - endif - open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1') - do a = 1, nV - do i = 1, nO - read(11,'(F20.12)') t1(i,a) + if (read_amplitudes) then + iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','r') + do a = 1, nV + do i = 1, nO + read(iunit,'(F20.12)') t1(i,a) + enddo enddo - enddo - close(11) + close(iunit) + endif end @@ -188,26 +186,23 @@ subroutine read_t2(nO,nV,t2) double precision, intent(out) :: t2(nO, nO, nV, nV) ! internal - integer :: i,j,a,b + integer :: i,j,a,b, iunit logical :: ok - inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok) - if (.not. ok) then - print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1' - print*, 'Do a first calculation with cc_write_t2 = True' - print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read' - call abort - endif - open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2') - do b = 1, nV - do a = 1, nV - do j = 1, nO - do i = 1, nO - read(11,'(F20.12)') t2(i,j,a,b) + integer, external :: getunitandopen + + if (read_amplitudes) then + iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','r') + do b = 1, nV + do a = 1, nV + do j = 1, nO + do i = 1, nO + read(iunit,'(F20.12)') t2(i,j,a,b) + enddo enddo enddo enddo - enddo - close(11) + close(iunit) + endif end From 2f246780eb0a08ee87ab5d98fe9c0e2f17685594 Mon Sep 17 00:00:00 2001 From: ydamour Date: Tue, 13 Jun 2023 14:05:13 +0200 Subject: [PATCH 3/4] fix bug in get_excitation_general --- src/utils_cc/org/phase.org | 2 ++ src/utils_cc/phase.irp.f | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/utils_cc/org/phase.org b/src/utils_cc/org/phase.org index 5f67859c..2156a251 100644 --- a/src/utils_cc/org/phase.org +++ b/src/utils_cc/org/phase.org @@ -137,6 +137,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n1(j) + if (k > n_anni(j)) exit if (l1(i,j) /= list_anni(k,j)) cycle pos_anni(k,j) = i k = k + 1 @@ -147,6 +148,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n2(j) + if (k > n_crea(j)) exit if (l2(i,j) /= list_crea(k,j)) cycle pos_crea(k,j) = i k = k + 1 diff --git a/src/utils_cc/phase.irp.f b/src/utils_cc/phase.irp.f index 01b41f49..e0703fb8 100644 --- a/src/utils_cc/phase.irp.f +++ b/src/utils_cc/phase.irp.f @@ -96,6 +96,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n1(j) + if (k > n_anni(j)) exit if (l1(i,j) /= list_anni(k,j)) cycle pos_anni(k,j) = i k = k + 1 @@ -106,6 +107,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n2(j) + if (k > n_crea(j)) exit if (l2(i,j) /= list_crea(k,j)) cycle pos_crea(k,j) = i k = k + 1 From a56644a808e0aea4c16d72c61b717ac4b2e0cabc Mon Sep 17 00:00:00 2001 From: ydamour Date: Tue, 13 Jun 2023 14:24:39 +0200 Subject: [PATCH 4/4] remove old stuffs --- src/mo_optimization/my_providers.irp.f | 141 ------------------------- 1 file changed, 141 deletions(-) delete mode 100644 src/mo_optimization/my_providers.irp.f diff --git a/src/mo_optimization/my_providers.irp.f b/src/mo_optimization/my_providers.irp.f deleted file mode 100644 index 7469ffd5..00000000 --- a/src/mo_optimization/my_providers.irp.f +++ /dev/null @@ -1,141 +0,0 @@ -! Dimensions of MOs - - -BEGIN_PROVIDER [ integer, n_mo_dim ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of MOs we can build, - ! with i>j - END_DOC - - n_mo_dim = mo_num*(mo_num-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_core ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of core MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_core = dim_list_core_orb*(dim_list_core_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_act ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of active MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_act = dim_list_act_orb*(dim_list_act_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_inact ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of inactive MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_inact = dim_list_inact_orb*(dim_list_inact_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_virt ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of virtual MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_virt = dim_list_virt_orb*(dim_list_virt_orb-1)/2 - -END_PROVIDER - -! Energies/criterions - -BEGIN_PROVIDER [ double precision, my_st_av_energy ] - implicit none - BEGIN_DOC - ! State average CI energy - END_DOC - - !call update_st_av_ci_energy(my_st_av_energy) - call state_average_energy(my_st_av_energy) - -END_PROVIDER - -! With all the MOs - -BEGIN_PROVIDER [ double precision, my_gradient_opt, (n_mo_dim) ] -&BEGIN_PROVIDER [ double precision, my_CC1_opt ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, for all the MOs. - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision :: norm_grad - - PROVIDE mo_two_e_integrals_in_map - - call gradient_opt(n_mo_dim, my_gradient_opt, my_CC1_opt, norm_grad) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, my_hessian_opt, (n_mo_dim, n_mo_dim) ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, for all the MOs. - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision, allocatable :: h_f(:,:,:,:) - - PROVIDE mo_two_e_integrals_in_map - - allocate(h_f(mo_num, mo_num, mo_num, mo_num)) - - call hessian_list_opt(n_mo_dim, my_hessian_opt, h_f) - -END_PROVIDER - -! With the list of active MOs -! Can be generalized to any mo_class by changing the list/dimension - -BEGIN_PROVIDER [ double precision, my_gradient_list_opt, (n_mo_dim_act) ] -&BEGIN_PROVIDER [ double precision, my_CC2_opt ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, only for the active MOs ! - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision :: norm_grad - - PROVIDE mo_two_e_integrals_in_map !one_e_dm_mo two_e_dm_mo mo_one_e_integrals - - call gradient_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_list_opt, my_CC2_opt, norm_grad) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, my_hessian_list_opt, (n_mo_dim_act, n_mo_dim_act) ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, only for the active MOs ! - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision, allocatable :: h_f(:,:,:,:) - - PROVIDE mo_two_e_integrals_in_map - - allocate(h_f(dim_list_act_orb, dim_list_act_orb, dim_list_act_orb, dim_list_act_orb)) - - call hessian_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_list_opt, h_f) - -END_PROVIDER