From 4c0de615fb0a572212c2c3940c601f64e8cd0164 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Apr 2023 19:36:39 +0200 Subject: [PATCH 01/20] Fix qp_extract_cipsi_data.py --- scripts/qp_extract_cipsi_data.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/qp_extract_cipsi_data.py b/scripts/qp_extract_cipsi_data.py index 8f0b1f3c..dd8e9c4d 100755 --- a/scripts/qp_extract_cipsi_data.py +++ b/scripts/qp_extract_cipsi_data.py @@ -23,6 +23,9 @@ def extract_data(output): reading = False for iline, line in enumerate(lines): + if line.startswith("Summary at N_det"): + reading = False + if not reading and line.startswith(" N_det "): n_det = int(re.search(r"N_det\s+=\s+(\d+)", line).group(1)) reading = True From 44d867297423cf33736340b111dc28f6f114a5dc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 17 Apr 2023 16:24:07 +0200 Subject: [PATCH 02/20] OMP Critial around format_w_error --- src/utils/format_w_error.irp.f | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/utils/format_w_error.irp.f b/src/utils/format_w_error.irp.f index 1378d367..ce2665a7 100644 --- a/src/utils/format_w_error.irp.f +++ b/src/utils/format_w_error.irp.f @@ -1,7 +1,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error) implicit none - + BEGIN_DOC ! Format for double precision, value(error) END_DOC @@ -14,7 +14,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err ! out ! | format_value | character | string FX.Y for the format | - ! | str_error | character | string of the error | + ! | str_error | character | string of the error | ! internal ! | str_size | character | size in string format | @@ -33,6 +33,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err character(len=20) :: str_size, str_nb_digits, str_exp integer :: nb_digits + !$OMP CRITICAL ! max_nb_digit: Y max ! size_nb = Size of the double: X (FX.Y) write(str_size,'(I3)') size_nb @@ -40,17 +41,17 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err ! Error write(str_exp,'(1pE20.0)') error str_error = trim(adjustl(str_exp)) - + ! Number of digit: Y (FX.Y) from the exponent str_nb_digits = str_exp(19:20) read(str_nb_digits,*) nb_digits - + ! If the error is 0d0 - if (error <= 1d-16) then + if (error <= 1d-16) then write(str_nb_digits,*) max_nb_digits endif - ! If the error is too small + ! If the error is too small if (nb_digits > max_nb_digits) then write(str_nb_digits,*) max_nb_digits str_error(1:1) = '0' @@ -65,7 +66,8 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err ! FX.Y,A1,A1,A1 for value(str_error) !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1' - ! FX.Y just for the value + ! FX.Y just for the value format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits)) + !$OMP END CRITICAL end From 5b6ecfa564b8d889981342c2e9ad597d124d15a3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 17 Apr 2023 17:03:16 +0200 Subject: [PATCH 03/20] Improve thread-safety --- src/cipsi/environment.irp.f | 2 ++ src/davidson/davidson_parallel.irp.f | 4 ++++ src/ezfio_files/lock.irp.f | 17 +++++++++++++++++ src/mo_two_e_ints/map_integrals.irp.f | 2 ++ src/nuclei/nuclei.irp.f | 14 +++++++------- src/two_body_rdm/io_two_rdm.irp.f | 12 ++++++++---- src/utils/format_w_error.irp.f | 4 ++-- src/utils/memory.irp.f | 8 ++++++-- 8 files changed, 48 insertions(+), 15 deletions(-) diff --git a/src/cipsi/environment.irp.f b/src/cipsi/environment.irp.f index 5c0e0820..363b8f1c 100644 --- a/src/cipsi/environment.irp.f +++ b/src/cipsi/environment.irp.f @@ -7,7 +7,9 @@ BEGIN_PROVIDER [ integer, nthreads_pt2 ] character*(32) :: env call getenv('QP_NTHREADS_PT2',env) if (trim(env) /= '') then + call lock_io() read(env,*) nthreads_pt2 + call unlock_io() call write_int(6,nthreads_pt2,'Target number of threads for PT2') endif END_PROVIDER diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index e627dfc9..399ab11b 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -150,7 +150,9 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze, exit endif if(task_id == 0) exit + call lock_io() read (msg,*) imin, imax, ishift, istep + call unlock_io() integer :: k do k=imin,imax v_t(:,k) = 0.d0 @@ -555,7 +557,9 @@ BEGIN_PROVIDER [ integer, nthreads_davidson ] character*(32) :: env call getenv('QP_NTHREADS_DAVIDSON',env) if (trim(env) /= '') then + call lock_io() read(env,*) nthreads_davidson + call unlock_io() call write_int(6,nthreads_davidson,'Target number of threads for ') endif END_PROVIDER diff --git a/src/ezfio_files/lock.irp.f b/src/ezfio_files/lock.irp.f index 53a99254..d28f7641 100644 --- a/src/ezfio_files/lock.irp.f +++ b/src/ezfio_files/lock.irp.f @@ -9,4 +9,21 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), file_lock ] call omp_init_lock(file_lock) END_PROVIDER +! These functions need to be called because internal read and write are not thread safe. +subroutine lock_io() + implicit none + BEGIN_DOC +! Needs to be called because before doing I/O because internal read and write +! are not thread safe. + END_DOC + call omp_set_lock(file_lock) +end subroutine lock_io() +subroutine unlock_io() + implicit none + BEGIN_DOC +! Needs to be called because afterdoing I/O because internal read and write +! are not thread safe. + END_DOC + call omp_unset_lock(file_lock) +end subroutine lock_io() diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 272916e3..ada256a2 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -377,6 +377,7 @@ integer function load_mo_integrals(filename) integer*8 :: n, j load_mo_integrals = 1 open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') + call lock_io() read(66,err=98,end=98) iknd, kknd if (iknd /= integral_kind) then print *, 'Wrong integrals kind in file :', iknd @@ -399,6 +400,7 @@ integer function load_mo_integrals(filename) n = mo_integrals_map%map(i)%n_elements read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) enddo + call unlock_io() call map_sort(mo_integrals_map) load_mo_integrals = 0 return diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index 3c04316f..fabdc42e 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -241,13 +241,13 @@ END_PROVIDER enddo character*(80) :: buffer, dummy do - read(iunit,'(A80)',end=10) buffer - read(buffer,*) i ! First read i - read(buffer,*) i, element_name(i), dummy, element_mass(i) - enddo - 10 continue - close(10) - endif + read(iunit,'(A80)',end=10) buffer + read(buffer,*) i ! First read i + read(buffer,*) i, element_name(i), dummy, element_mass(i) + enddo + 10 continue + close(10) + endif IRP_IF MPI_DEBUG print *, irp_here, mpi_rank diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f index f7008ca9..bdd8a4f9 100644 --- a/src/two_body_rdm/io_two_rdm.irp.f +++ b/src/two_body_rdm/io_two_rdm.irp.f @@ -1,15 +1,17 @@ subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file) implicit none integer, intent(in) :: n_orb,nstates - character*(128), intent(in) :: name_file + character*(128), intent(in) :: name_file double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates) character*(128) :: output integer :: i_unit_output,getUnitAndOpen - PROVIDE ezfio_filename + PROVIDE ezfio_filename output=trim(ezfio_filename)//'/work/'//trim(name_file) i_unit_output = getUnitAndOpen(output,'W') + call lock_io() write(i_unit_output)array_tmp + call unlock_io() close(unit=i_unit_output) end @@ -18,12 +20,14 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file) character*(128) :: output integer :: i_unit_output,getUnitAndOpen integer, intent(in) :: n_orb,nstates - character*(128), intent(in) :: name_file + character*(128), intent(in) :: name_file double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states) - PROVIDE ezfio_filename + PROVIDE ezfio_filename output=trim(ezfio_filename)//'/work/'//trim(name_file) i_unit_output = getUnitAndOpen(output,'R') + call lock_io() read(i_unit_output)array_tmp + call unlock_io() close(unit=i_unit_output) end diff --git a/src/utils/format_w_error.irp.f b/src/utils/format_w_error.irp.f index ce2665a7..7f7458b6 100644 --- a/src/utils/format_w_error.irp.f +++ b/src/utils/format_w_error.irp.f @@ -33,7 +33,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err character(len=20) :: str_size, str_nb_digits, str_exp integer :: nb_digits - !$OMP CRITICAL + call lock_io() ! max_nb_digit: Y max ! size_nb = Size of the double: X (FX.Y) write(str_size,'(I3)') size_nb @@ -68,6 +68,6 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err ! FX.Y just for the value format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits)) - !$OMP END CRITICAL + call unlock_io() end diff --git a/src/utils/memory.irp.f b/src/utils/memory.irp.f index d5a066a1..115b2cbe 100644 --- a/src/utils/memory.irp.f +++ b/src/utils/memory.irp.f @@ -8,7 +8,9 @@ BEGIN_PROVIDER [ integer, qp_max_mem ] qp_max_mem = 2000 call getenv('QP_MAXMEM',env) if (trim(env) /= '') then + call lock_io() read(env,*) qp_max_mem + call unlock_io() endif call write_int(6,qp_max_mem,'Target maximum memory (GB)') @@ -25,7 +27,7 @@ subroutine resident_memory(value) character*(32) :: key double precision, intent(out) :: value - call omp_set_lock(file_lock) + call lock_io() call usleep(10) value = 0.d0 @@ -40,7 +42,7 @@ subroutine resident_memory(value) 20 continue close(iunit) value = value / (1024.d0*1024.d0) - call omp_unset_lock(file_lock) + call unlock_io() end function subroutine total_memory(value) @@ -53,6 +55,7 @@ subroutine total_memory(value) character*(32) :: key double precision, intent(out) :: value + call lock_io() iunit = getUnitAndOpen('/proc/self/status','r') do read(iunit,*,err=10,end=20) key, value @@ -64,6 +67,7 @@ subroutine total_memory(value) 20 continue close(iunit) value = value / (1024.d0*1024.d0) + call unlock_io() end function double precision function memory_of_double(n) From 6b4bf5b601fb759772be99af301b4959d17c18e6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 18 Apr 2023 14:41:08 +0200 Subject: [PATCH 04/20] Raise error --- external/ezfio | 2 +- external/irpf90 | 2 +- external/qp2-dependencies | 2 +- scripts/compilation/qp_create_ninja | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/external/ezfio b/external/ezfio index d5805497..ed1df9f3 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 +Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c diff --git a/external/irpf90 b/external/irpf90 index 0007f72f..33ca5e10 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 +Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 6e23ebac..9e5b27ce 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 6e23ebac001acae91d1c762ca934e09a9b7d614a +Subproject commit 9e5b27ce5a174901765cec9db9e7b2aa6170a5de diff --git a/scripts/compilation/qp_create_ninja b/scripts/compilation/qp_create_ninja index aad85778..27b34901 100755 --- a/scripts/compilation/qp_create_ninja +++ b/scripts/compilation/qp_create_ninja @@ -25,7 +25,7 @@ except ImportError: "quantum_package.rc")) print("\n".join(["", "Error:", "source %s" % f, ""])) - sys.exit(1) + raise # Compress path def comp_path(path): From 750ec9ca00559eced5f08ece573bd583f5aa8c2b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 19 Apr 2023 09:31:36 +0200 Subject: [PATCH 05/20] Fix f77zmq --- configure | 2 +- external/qp2-dependencies | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/configure b/configure index 5b50d0d7..d3377093 100755 --- a/configure +++ b/configure @@ -232,7 +232,7 @@ EOF execute << EOF cd "\${QP_ROOT}"/external - tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.2.tar.gz + tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.?.tar.gz cd f77-zmq-* ./configure --prefix=\$QP_ROOT export ZMQ_H="\$QP_ROOT"/include/zmq.h diff --git a/external/qp2-dependencies b/external/qp2-dependencies index ce14f57b..fd43778e 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 +Subproject commit fd43778e12bb5858c4c780c34346be0f158b8cc7 From 94662d3da070a8f4602578a62b3734fcce709593 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 Apr 2023 15:11:50 +0200 Subject: [PATCH 06/20] Introduced JSON --- src/hartree_fock/NEED | 1 + src/hartree_fock/scf.irp.f | 6 +++++ src/json/EZFIO.cfg | 5 ++++ src/json/NEED | 1 + src/json/README.rst | 5 ++++ src/json/json.irp.f | 39 +++++++++++++++++++++++++++ src/json/json_formats.irp.f | 26 ++++++++++++++++++ src/scf_utils/roothaan_hall_scf.irp.f | 39 ++++++++++++++++++++++++--- src/two_body_rdm/io_two_rdm.irp.f | 8 +++--- 9 files changed, 122 insertions(+), 8 deletions(-) create mode 100644 src/json/EZFIO.cfg create mode 100644 src/json/NEED create mode 100644 src/json/README.rst create mode 100644 src/json/json.irp.f create mode 100644 src/json/json_formats.irp.f diff --git a/src/hartree_fock/NEED b/src/hartree_fock/NEED index 2b3fa238..e168bd80 100644 --- a/src/hartree_fock/NEED +++ b/src/hartree_fock/NEED @@ -1,3 +1,4 @@ ao_one_e_ints ao_two_e_ints scf_utils +json diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index a7ac9fe4..d22b11ab 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -80,9 +80,15 @@ subroutine run mo_label = 'Orthonormalized' + write(json_unit,*) '"scf" : [' + call Roothaan_Hall_SCF call ezfio_set_hartree_fock_energy(SCF_energy) + write(json_unit,*) ']' + + call json_close + end diff --git a/src/json/EZFIO.cfg b/src/json/EZFIO.cfg new file mode 100644 index 00000000..7dc8d796 --- /dev/null +++ b/src/json/EZFIO.cfg @@ -0,0 +1,5 @@ +[empty] +type: logical +doc: Needed to create the json directory +interface: ezfio + diff --git a/src/json/NEED b/src/json/NEED new file mode 100644 index 00000000..5a3182ed --- /dev/null +++ b/src/json/NEED @@ -0,0 +1 @@ +ezfio_files diff --git a/src/json/README.rst b/src/json/README.rst new file mode 100644 index 00000000..3dd9ffbb --- /dev/null +++ b/src/json/README.rst @@ -0,0 +1,5 @@ +==== +json +==== + +JSON files to simplify getting output information from QP. diff --git a/src/json/json.irp.f b/src/json/json.irp.f new file mode 100644 index 00000000..5a92f22f --- /dev/null +++ b/src/json/json.irp.f @@ -0,0 +1,39 @@ +BEGIN_PROVIDER [ character*(128), json_filename ] + implicit none + BEGIN_DOC + ! Fortran unit of the JSON file + END_DOC + integer, external :: getUnitAndOpen + integer :: counter + character*(128) :: prefix + logical :: exists + + prefix = trim(ezfio_filename)//'/json/' + + exists = .True. + counter = 0 + do while (exists) + counter += 1 + write(json_filename, '(A,I5.5,A)') trim(prefix), counter, '.json' + INQUIRE(FILE=trim(json_filename), EXIST=exists) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ integer, json_unit] + implicit none + BEGIN_DOC + ! Unit file for JSON output + END_DOC + integer, external :: getUnitAndOpen + call ezfio_set_json_empty(.False.) + json_unit = getUnitAndOpen(json_filename, 'w') + write(json_unit, '(A)') '{' +END_PROVIDER + +subroutine json_close + write(json_unit, '(A)') '}' + close(json_unit) + FREE json_unit +end + diff --git a/src/json/json_formats.irp.f b/src/json/json_formats.irp.f new file mode 100644 index 00000000..14a8f014 --- /dev/null +++ b/src/json/json_formats.irp.f @@ -0,0 +1,26 @@ + BEGIN_PROVIDER [ character*(64), json_int_fmt ] +&BEGIN_PROVIDER [ character*(64), json_int_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_real_fmt ] +&BEGIN_PROVIDER [ character*(64), json_real_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_str_fmt ] +&BEGIN_PROVIDER [ character*(64), json_str_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_true_fmt ] +&BEGIN_PROVIDER [ character*(64), json_true_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_false_fmt ] +&BEGIN_PROVIDER [ character*(64), json_false_fmtx ] + implicit none + BEGIN_DOC + ! Formats for JSON output. + ! x: used to mark the last write (no comma) + END_DOC + json_int_fmt = '('' "'',A,''": '',I10,'','')' + json_int_fmtx = '('' "'',A,''": '',I10)' + json_real_fmt = '('' "'',A,''": '',E22.15,'','')' + json_real_fmtx = '('' "'',A,''": '',E22.15)' + json_str_fmt = '('' "'',A,''": "'',A,''",'')' + json_str_fmtx = '('' "'',A,''": "'',A,''"'')' + json_true_fmt = '('' "'',A,''": true,'')' + json_true_fmtx = '('' "'',A,''": true'')' + json_false_fmt = '('' "'',A,''": false,'')' + json_false_fmtx = '('' "'',A,''": false'')' +END_PROVIDER diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 3b9eaeb4..449afdc8 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -12,6 +12,7 @@ END_DOC integer :: iteration_SCF,dim_DIIS,index_dim_DIIS + logical :: converged integer :: i,j logical, external :: qp_stop double precision, allocatable :: mo_coef_save(:,:) @@ -50,10 +51,8 @@ END_DOC ! PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO - do while ( & - ( (max_error_DIIS > threshold_DIIS_nonzero) .or. & - (dabs(Delta_energy_SCF) > thresh_SCF) & - ) .and. (iteration_SCF < n_it_SCF_max) ) + converged = .False. + do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) ) ! Increment cycle number @@ -144,17 +143,45 @@ END_DOC SOFT_TOUCH level_shift energy_SCF_previous = energy_SCF + converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. & + (dabs(Delta_energy_SCF) <= thresh_SCF) ) + ! Print results at the end of each iteration write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') & iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS +! Write data in JSON file + + if (iteration_SCF == 1) then + write(json_unit, *) '{' + else + write(json_unit, *) '}, {' + endif + write(json_unit, json_int_fmt) 'iteration', iteration_SCF + write(json_unit, json_real_fmt) 'energy', energy_SCF + write(json_unit, json_real_fmt) 'delta_energy_SCF', Delta_energy_SCF + write(json_unit, json_real_fmt) 'max_error_DIIS', max_error_DIIS + write(json_unit, json_real_fmt) 'level_shift', level_shift + write(json_unit, json_int_fmt) 'dim_DIIS', dim_DIIS + if (Delta_energy_SCF < 0.d0) then call save_mos + write(json_unit, json_true_fmt) 'saved' + else + write(json_unit, json_false_fmt) 'saved' endif + + if (converged) then + write(json_unit, json_true_fmtx) 'converged' + else + write(json_unit, json_false_fmtx) 'converged' + endif + if (qp_stop()) exit enddo + write(json_unit, *) '}' if (iteration_SCF < n_it_SCF_max) then mo_label = 'Canonical' @@ -166,6 +193,10 @@ END_DOC write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & '====','================','================','================','================' write(6,*) + if (converged) then + write(6,*) 'SCF converged' + endif + if(.not.frozen_orb_scf)then call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), & diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f index f7008ca9..67837682 100644 --- a/src/two_body_rdm/io_two_rdm.irp.f +++ b/src/two_body_rdm/io_two_rdm.irp.f @@ -1,12 +1,12 @@ subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file) implicit none integer, intent(in) :: n_orb,nstates - character*(128), intent(in) :: name_file + character*(128), intent(in) :: name_file double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates) character*(128) :: output integer :: i_unit_output,getUnitAndOpen - PROVIDE ezfio_filename + PROVIDE ezfio_filename output=trim(ezfio_filename)//'/work/'//trim(name_file) i_unit_output = getUnitAndOpen(output,'W') write(i_unit_output)array_tmp @@ -18,9 +18,9 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file) character*(128) :: output integer :: i_unit_output,getUnitAndOpen integer, intent(in) :: n_orb,nstates - character*(128), intent(in) :: name_file + character*(128), intent(in) :: name_file double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states) - PROVIDE ezfio_filename + PROVIDE ezfio_filename output=trim(ezfio_filename)//'/work/'//trim(name_file) i_unit_output = getUnitAndOpen(output,'R') read(i_unit_output)array_tmp From 5039bb674d6d2433de182011bbb8973748588c25 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 Apr 2023 18:06:37 +0200 Subject: [PATCH 07/20] Fixed need for JSON --- external/ezfio | 2 +- external/irpf90 | 2 +- external/qp2-dependencies | 2 +- src/json/json.irp.f | 6 ++++++ src/kohn_sham/ks_scf.irp.f | 4 ++++ src/kohn_sham_rs/rs_ks_scf.irp.f | 3 +++ src/scf_utils/NEED | 1 + src/scf_utils/roothaan_hall_scf.irp.f | 4 ++++ 8 files changed, 21 insertions(+), 3 deletions(-) diff --git a/external/ezfio b/external/ezfio index ed1df9f3..d5805497 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c +Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 diff --git a/external/irpf90 b/external/irpf90 index 33ca5e10..0007f72f 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 +Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 diff --git a/external/qp2-dependencies b/external/qp2-dependencies index fd43778e..e0d0e02e 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit fd43778e12bb5858c4c780c34346be0f158b8cc7 +Subproject commit e0d0e02e9f5ece138d1520106954a881ab0b8db2 diff --git a/src/json/json.irp.f b/src/json/json.irp.f index 5a92f22f..1fc24eb6 100644 --- a/src/json/json.irp.f +++ b/src/json/json.irp.f @@ -10,6 +10,7 @@ BEGIN_PROVIDER [ character*(128), json_filename ] prefix = trim(ezfio_filename)//'/json/' + call lock_io exists = .True. counter = 0 do while (exists) @@ -17,6 +18,7 @@ BEGIN_PROVIDER [ character*(128), json_filename ] write(json_filename, '(A,I5.5,A)') trim(prefix), counter, '.json' INQUIRE(FILE=trim(json_filename), EXIST=exists) enddo + call unlock_io END_PROVIDER @@ -27,13 +29,17 @@ BEGIN_PROVIDER [ integer, json_unit] END_DOC integer, external :: getUnitAndOpen call ezfio_set_json_empty(.False.) + call lock_io json_unit = getUnitAndOpen(json_filename, 'w') write(json_unit, '(A)') '{' + call unlock_io END_PROVIDER subroutine json_close + call lock_io write(json_unit, '(A)') '}' close(json_unit) + call unlock_io FREE json_unit end diff --git a/src/kohn_sham/ks_scf.irp.f b/src/kohn_sham/ks_scf.irp.f index aa6efd52..85bfc333 100644 --- a/src/kohn_sham/ks_scf.irp.f +++ b/src/kohn_sham/ks_scf.irp.f @@ -90,7 +90,11 @@ subroutine run ! Choose SCF algorithm + write(json_unit,*) '"scf" : [' call Roothaan_Hall_SCF + write(json_unit,*) ']' + + call json_close end diff --git a/src/kohn_sham_rs/rs_ks_scf.irp.f b/src/kohn_sham_rs/rs_ks_scf.irp.f index 84b85136..f28fd861 100644 --- a/src/kohn_sham_rs/rs_ks_scf.irp.f +++ b/src/kohn_sham_rs/rs_ks_scf.irp.f @@ -93,7 +93,10 @@ subroutine run level_shift += 1.d0 touch level_shift + write(json_unit,*) '"scf" : [' call Roothaan_Hall_SCF + write(json_unit,*) ']' + call json_close call ezfio_set_kohn_sham_rs_energy(SCF_energy) write(*, '(A22,X,F16.10)') 'one_e_energy = ',one_e_energy diff --git a/src/scf_utils/NEED b/src/scf_utils/NEED index b89695da..292d343a 100644 --- a/src/scf_utils/NEED +++ b/src/scf_utils/NEED @@ -1,2 +1,3 @@ mo_guess bitmask +json diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 449afdc8..08fe7acf 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -153,6 +153,7 @@ END_DOC ! Write data in JSON file + call lock_io if (iteration_SCF == 1) then write(json_unit, *) '{' else @@ -164,6 +165,7 @@ END_DOC write(json_unit, json_real_fmt) 'max_error_DIIS', max_error_DIIS write(json_unit, json_real_fmt) 'level_shift', level_shift write(json_unit, json_int_fmt) 'dim_DIIS', dim_DIIS + call unlock_io if (Delta_energy_SCF < 0.d0) then call save_mos @@ -172,11 +174,13 @@ END_DOC write(json_unit, json_false_fmt) 'saved' endif + call lock_io if (converged) then write(json_unit, json_true_fmtx) 'converged' else write(json_unit, json_false_fmtx) 'converged' endif + call unlock_io if (qp_stop()) exit From 528bf20e1e00e2677f695bba8b61c203cc053777 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 Apr 2023 18:10:08 +0200 Subject: [PATCH 08/20] Cleaning useless function --- src/kohn_sham/print_mos.irp.f | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/kohn_sham/print_mos.irp.f b/src/kohn_sham/print_mos.irp.f index 5e728444..7105c989 100644 --- a/src/kohn_sham/print_mos.irp.f +++ b/src/kohn_sham/print_mos.irp.f @@ -21,10 +21,3 @@ program print_mos end -double precision function f_mu(x) - implicit none - double precision, intent(in) :: x - - - -end From 5fb6ed0180090e1d7d5cd009ce0e3588815670a3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 00:50:07 +0200 Subject: [PATCH 09/20] Added JSON in FCI --- src/cipsi/NEED | 1 + src/cipsi/cipsi.irp.f | 10 +-- src/cipsi/stochastic_cipsi.irp.f | 8 +-- src/cipsi/write_cipsi_json.irp.f | 53 ++++++++++++++ src/fci/fci.irp.f | 7 ++ src/hartree_fock/scf.irp.f | 7 +- src/iterations/EZFIO.cfg | 24 ------- src/iterations/io.irp.f | 37 ---------- src/iterations/iterations.irp.f | 91 +++++++++++++++--------- src/iterations/print_extrapolation.irp.f | 14 ++-- src/json/json_formats.irp.f | 20 ++++++ src/scf_utils/roothaan_hall_scf.irp.f | 6 +- 12 files changed, 162 insertions(+), 116 deletions(-) create mode 100644 src/cipsi/write_cipsi_json.irp.f delete mode 100644 src/iterations/EZFIO.cfg delete mode 100644 src/iterations/io.irp.f diff --git a/src/cipsi/NEED b/src/cipsi/NEED index 5bd742bc..89c128ec 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -1,3 +1,4 @@ +json perturbation zmq mpi diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 5225c6df..88aaeae0 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -16,7 +16,6 @@ subroutine run_cipsi double precision, external :: memory_of_double PROVIDE H_apply_buffer_allocated - N_iter = 1 threshold_generators = 1.d0 SOFT_TOUCH threshold_generators @@ -76,7 +75,6 @@ subroutine run_cipsi ) write(*,'(A)') '--------------------------------------------------------------------------------' - to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) if (do_pt2) then @@ -106,10 +104,10 @@ subroutine run_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) + call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) call print_extrapolated_energy() call print_mol_properties() - N_iter += 1 + call write_cipsi_json(pt2_data,pt2_data_err) if (qp_stop()) exit @@ -155,11 +153,13 @@ subroutine run_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep(1:N_states), & pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) + call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) call print_extrapolated_energy() call print_mol_properties() + call write_cipsi_json(pt2_data,pt2_data_err) endif call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) end + diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 35e80eb8..b83e658a 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -15,7 +15,6 @@ subroutine run_stochastic_cipsi double precision, external :: memory_of_double PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map - N_iter = 1 threshold_generators = 1.d0 SOFT_TOUCH threshold_generators @@ -96,10 +95,10 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) + call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) call print_extrapolated_energy() call print_mol_properties() - N_iter += 1 + call write_cipsi_json(pt2_data,pt2_data_err) if (qp_stop()) exit @@ -135,9 +134,10 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep, & pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) + call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) call print_extrapolated_energy() call print_mol_properties() + call write_cipsi_json(pt2_data,pt2_data_err) endif call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/src/cipsi/write_cipsi_json.irp.f b/src/cipsi/write_cipsi_json.irp.f new file mode 100644 index 00000000..98a402a2 --- /dev/null +++ b/src/cipsi/write_cipsi_json.irp.f @@ -0,0 +1,53 @@ +subroutine write_cipsi_json(pt2_data, pt2_data_err) + use selection_types + implicit none + BEGIN_DOC +! Writes JSON data for CIPSI runs + END_DOC + type(pt2_type), intent(in) :: pt2_data, pt2_data_err + integer :: i,j,k + + call lock_io + character*(64), allocatable :: fmtk(:) + integer :: N_states_p, N_iter_p + N_states_p = min(N_states,N_det) + N_iter_p = min(N_iter,8) + allocate(fmtk(0:N_iter_p)) + fmtk(:) = '('' '',E22.15,'','')' + fmtk(N_iter_p) = '('' '',E22.15)' + + write(json_unit, json_dict_uopen_fmt) + write(json_unit, json_int_fmt) 'n_det', N_det + if (s2_eig) then + write(json_unit, json_int_fmt) 'n_cfg', N_configuration + if (only_expected_s2) then + write(json_unit, json_int_fmt) 'n_csf', N_csf + endif + endif + write(json_unit, json_array_open_fmt) 'states' + do k=1,N_states_p + write(json_unit, json_dict_uopen_fmt) + write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k) + write(json_unit, json_real_fmt) 's2', psi_s2(k) + write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k) + write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k) + write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k) + write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k) + write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k) + write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k) + write(json_unit, json_array_open_fmt) 'ex_energy' + do i=2,N_iter_p + write(json_unit, fmtk(i)) extrapolated_energy(i,k) + enddo + write(json_unit, json_array_close_fmtx) + if (k < N_states_p) then + write(json_unit, json_dict_close_fmt) + else + write(json_unit, json_dict_close_fmtx) + endif + enddo + write(json_unit, json_array_close_fmtx) + write(json_unit, json_dict_close_fmt) + deallocate(fmtk) + call unlock_io +end diff --git a/src/fci/fci.irp.f b/src/fci/fci.irp.f index 9d9c0b7d..bb2a93f8 100644 --- a/src/fci/fci.irp.f +++ b/src/fci/fci.irp.f @@ -39,12 +39,19 @@ program fci if (.not.is_zmq_slave) then PROVIDE psi_det psi_coef mo_two_e_integrals_in_map + write(json_unit,json_array_open_fmt) 'fci' + if (do_pt2) then call run_stochastic_cipsi else call run_cipsi endif + write(json_unit,json_dict_uopen_fmt) + write(json_unit,json_dict_close_fmtx) + write(json_unit,json_array_close_fmtx) + call json_close + else PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index d22b11ab..a361c04f 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -80,15 +80,14 @@ subroutine run mo_label = 'Orthonormalized' - write(json_unit,*) '"scf" : [' + write(json_unit,json_array_open_fmt) 'scf' call Roothaan_Hall_SCF - call ezfio_set_hartree_fock_energy(SCF_energy) - - write(json_unit,*) ']' + write(json_unit,json_array_close_fmtx) call json_close + call ezfio_set_hartree_fock_energy(SCF_energy) end diff --git a/src/iterations/EZFIO.cfg b/src/iterations/EZFIO.cfg deleted file mode 100644 index 2a5e94a7..00000000 --- a/src/iterations/EZFIO.cfg +++ /dev/null @@ -1,24 +0,0 @@ -[n_iter] -interface: ezfio -doc: Number of saved iterations -type:integer -default: 1 - -[n_det_iterations] -interface: ezfio, provider -doc: Number of determinants at each iteration -type: integer -size: (100) - -[energy_iterations] -interface: ezfio, provider -doc: The variational energy at each iteration -type: double precision -size: (determinants.n_states,100) - -[pt2_iterations] -interface: ezfio, provider -doc: The |PT2| correction at each iteration -type: double precision -size: (determinants.n_states,100) - diff --git a/src/iterations/io.irp.f b/src/iterations/io.irp.f deleted file mode 100644 index 821f5e84..00000000 --- a/src/iterations/io.irp.f +++ /dev/null @@ -1,37 +0,0 @@ -BEGIN_PROVIDER [ integer, n_iter ] - implicit none - BEGIN_DOC -! number of iterations - END_DOC - - logical :: has - PROVIDE ezfio_filename - if (mpi_master) then - - double precision :: zeros(N_states,100) - integer :: izeros(100) - zeros = 0.d0 - izeros = 0 - call ezfio_set_iterations_n_iter(0) - call ezfio_set_iterations_energy_iterations(zeros) - call ezfio_set_iterations_pt2_iterations(zeros) - call ezfio_set_iterations_n_det_iterations(izeros) - n_iter = 1 - endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read n_iter with MPI' - endif - IRP_ENDIF - - call write_time(6) - -END_PROVIDER - diff --git a/src/iterations/iterations.irp.f b/src/iterations/iterations.irp.f index 2c9faaf8..d06d1b6e 100644 --- a/src/iterations/iterations.irp.f +++ b/src/iterations/iterations.irp.f @@ -1,42 +1,65 @@ -BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ] - implicit none - BEGIN_DOC - ! Extrapolated energy, using E_var = f(PT2) where PT2=0 - END_DOC - integer :: i - do i=1,min(N_states,N_det) - call extrapolate_data(N_iter, & - energy_iterations(i,1:N_iter), & - pt2_iterations(i,1:N_iter), & - extrapolated_energy(1:N_iter,i)) - enddo -END_PROVIDER - - -subroutine save_iterations(e_, pt2_,n_) +BEGIN_PROVIDER [ integer, N_iter ] implicit none BEGIN_DOC -! Update the energy in the EZFIO file. +! Number of CIPSI iterations END_DOC - integer, intent(in) :: n_ - double precision, intent(in) :: e_(N_states), pt2_(N_states) - integer :: i - if (N_iter == 101) then - do i=2,N_iter-1 - energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter) - pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) + N_iter = 0 +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_iter_max ] + implicit none + BEGIN_DOC + ! Max number of iterations for extrapolations + END_DOC + N_iter_max = 8 +END_PROVIDER + + BEGIN_PROVIDER [ double precision, energy_iterations , (n_states,N_iter_max) ] +&BEGIN_PROVIDER [ double precision, pt2_iterations , (n_states,N_iter_max) ] +&BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter_max,N_states) ] + implicit none + BEGIN_DOC +! The energy at each iteration for the extrapolations + END_DOC + + energy_iterations = 0.d0 + pt2_iterations = 0.d0 + extrapolated_energy = 0.d0 +END_PROVIDER + +subroutine increment_n_iter(e, pt2_data) + use selection_types + implicit none + BEGIN_DOC +! Does what is necessary to increment n_iter + END_DOC + double precision, intent(in) :: e(*) + type(pt2_type), intent(in) :: pt2_data + integer :: k, i + + if (N_det < N_states) return + + if (N_iter < N_iter_max) then + N_iter += 1 + else + do k=2,N_iter + energy_iterations(1:N_states,k-1) = energy_iterations(1:N_states,k) + pt2_iterations(1:N_states,k-1) = pt2_iterations(1:N_states,k) enddo - N_iter = N_iter-1 - TOUCH N_iter endif + energy_iterations(1:N_states,N_iter) = e(1:N_states) + pt2_iterations(1:N_states,N_iter) = pt2_data % rpt2(1:N_states) - energy_iterations(1:N_states,N_iter) = e_(1:N_states) - pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states) - n_det_iterations(N_iter) = n_ - call ezfio_set_iterations_N_iter(N_iter) - call ezfio_set_iterations_energy_iterations(energy_iterations) - call ezfio_set_iterations_pt2_iterations(pt2_iterations) - call ezfio_set_iterations_n_det_iterations(n_det_iterations) + if (N_iter < 2) then + extrapolated_energy(1,:) = energy_iterations(:,1) + pt2_iterations(:,1) + extrapolated_energy(2,:) = energy_iterations(:,2) + pt2_iterations(:,2) + else + do i=1,N_states + call extrapolate_data(N_iter, & + energy_iterations(i,1:N_iter), & + pt2_iterations(i,1:N_iter), & + extrapolated_energy(1:N_iter,i)) + enddo + endif end - diff --git a/src/iterations/print_extrapolation.irp.f b/src/iterations/print_extrapolation.irp.f index cb46fb67..111429bf 100644 --- a/src/iterations/print_extrapolation.irp.f +++ b/src/iterations/print_extrapolation.irp.f @@ -5,10 +5,14 @@ subroutine print_extrapolated_energy END_DOC integer :: i,k + integer :: N_states_p, N_iter_p if (N_iter< 2) then return endif + N_states_p = min(N_states,N_det) + N_iter_p = min(N_iter, 8) + write(*,'(A)') '' write(*,'(A)') 'Extrapolated energies' write(*,'(A)') '------------------------' @@ -20,20 +24,20 @@ subroutine print_extrapolated_energy write(*,*) '=========== ', '===================' write(*,*) 'minimum PT2 ', 'Extrapolated energy' write(*,*) '=========== ', '===================' - do k=2,min(N_iter,8) - write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1) + do k=2,N_iter_p + write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,k), extrapolated_energy(k,1) enddo write(*,*) '=========== ', '===================' - do i=2, min(N_states,N_det) + do i=2, N_states_p print *, '' print *, 'State ', i print *, '' write(*,*) '=========== ', '=================== ', '=================== ', '===================' write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) ' write(*,*) '=========== ', '=================== ', '=================== ', '===================' - do k=2,min(N_iter,8) - write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), & + do k=2,N_iter_p + write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), & extrapolated_energy(k,i) - extrapolated_energy(k,1), & (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0 enddo diff --git a/src/json/json_formats.irp.f b/src/json/json_formats.irp.f index 14a8f014..773114ba 100644 --- a/src/json/json_formats.irp.f +++ b/src/json/json_formats.irp.f @@ -8,6 +8,16 @@ &BEGIN_PROVIDER [ character*(64), json_true_fmtx ] &BEGIN_PROVIDER [ character*(64), json_false_fmt ] &BEGIN_PROVIDER [ character*(64), json_false_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_array_open_fmt ] +&BEGIN_PROVIDER [ character*(64), json_array_uopen_fmt ] +&BEGIN_PROVIDER [ character*(64), json_array_close_fmt ] +&BEGIN_PROVIDER [ character*(64), json_array_close_uopen_fmt ] +&BEGIN_PROVIDER [ character*(64), json_array_close_fmtx ] +&BEGIN_PROVIDER [ character*(64), json_dict_open_fmt ] +&BEGIN_PROVIDER [ character*(64), json_dict_uopen_fmt ] +&BEGIN_PROVIDER [ character*(64), json_dict_close_uopen_fmt ] +&BEGIN_PROVIDER [ character*(64), json_dict_close_fmt ] +&BEGIN_PROVIDER [ character*(64), json_dict_close_fmtx ] implicit none BEGIN_DOC ! Formats for JSON output. @@ -23,4 +33,14 @@ json_true_fmtx = '('' "'',A,''": true'')' json_false_fmt = '('' "'',A,''": false,'')' json_false_fmtx = '('' "'',A,''": false'')' + json_array_open_fmt = '('' "'',A,''": ['')' + json_array_uopen_fmt = '('' ['')' + json_array_close_fmt = '('' ],'')' + json_array_close_uopen_fmt = '('' ], ['')' + json_array_close_fmtx = '('' ]'')' + json_dict_open_fmt = '('' "'',A,''": {'')' + json_dict_uopen_fmt = '('' {'')' + json_dict_close_fmt = '('' },'')' + json_dict_close_uopen_fmt = '('' }, {'')' + json_dict_close_fmtx = '('' }'')' END_PROVIDER diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 08fe7acf..cf006035 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -155,9 +155,9 @@ END_DOC call lock_io if (iteration_SCF == 1) then - write(json_unit, *) '{' + write(json_unit, json_dict_uopen_fmt) else - write(json_unit, *) '}, {' + write(json_unit, json_dict_close_uopen_fmt) endif write(json_unit, json_int_fmt) 'iteration', iteration_SCF write(json_unit, json_real_fmt) 'energy', energy_SCF @@ -185,7 +185,7 @@ END_DOC if (qp_stop()) exit enddo - write(json_unit, *) '}' + write(json_unit, json_dict_close_fmtx) if (iteration_SCF < n_it_SCF_max) then mo_label = 'Canonical' From dd5291d90dca05f50976b10d42d9da5657641058 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 01:01:31 +0200 Subject: [PATCH 10/20] Added exc_energy_error.py --- scripts/exc_energy_error.py | 186 ++++++++++++++++++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100755 scripts/exc_energy_error.py diff --git a/scripts/exc_energy_error.py b/scripts/exc_energy_error.py new file mode 100755 index 00000000..ba9d7917 --- /dev/null +++ b/scripts/exc_energy_error.py @@ -0,0 +1,186 @@ +#!/usr/bin/env python +# Computes the error on the excitation energy of a CIPSI run. + +def student(p,df): + import scipy + from scipy.stats import t + return t.ppf(p, df) + + +def chi2cdf(x,k): + import scipy + import scipy.stats + return scipy.stats.chi2.cdf(x,k) + + +def jarque_bera(data): + + n = max(len(data), 2) + norm = 1./ sum( [ w for (_,w) in data ] ) + + mu = sum( [ w* x for (x,w) in data ] ) * norm + sigma2 = sum( [ w*(x-mu)**2 for (x,w) in data ] ) * norm + if sigma2 > 0.: + S = ( sum( [ w*(x-mu)**3 for (x,w) in data ] ) * norm ) / sigma2**(3./2.) + K = ( sum( [ w*(x-mu)**4 for (x,w) in data ] ) * norm ) / sigma2**2 + else: + S = 0. + K = 0. + + # Value of the Jarque-Bera test + JB = n/6. * (S**2 + 1./4. * (K-3.)**2) + + # Probability that the data comes from a Gaussian distribution + p = 1. - chi2cdf(JB,2) + + return JB, mu, sqrt(sigma2/(n-1)), p + + + +to_eV = 27.2107362681 +import sys, os +import scipy +import scipy.stats +from math import sqrt, gamma, exp +import json + + +def read_data(filename,state): + """ Read energies and PT2 from input file """ + with open(filename,'r') as f: + lines = json.load(f)['fci'] + + print(f"State: {state}") + + gs = [] + es = [] + + for l in lines: + try: + pt2_0 = l['states'][0]['pt2'] + e_0 = l['states'][0]['energy'] + pt2_1 = l['states'][state]['pt2'] + e_1 = l['states'][state]['energy'] + gs.append( (e_0, pt2_0) ) + es.append( (e_1, pt2_1) ) + except: pass + + def f(p_1, p0, p1): + e, pt2 = p0 + y0, x0 = p_1 + y1, x1 = p1 + try: + alpha = (y1-y0)/(x0-x1) + except ZeroDivisionError: + alpha = 1. + return [e, pt2, alpha] + + for l in (gs, es): + p_1, p0, p1 = l[0], l[0], l[1] + l[0] = f(p_1, p0, p1) + + for i in range(1,len(l)-1): + p_1 = (l[i-1][0], l[i-1][1]) + p0 = l[i] + p1 = l[i+1] + l[i] = f(p_1, p0, p1) + + i = len(l)-1 + p_1 = (l[i-1][0], l[i-1][1]) + p0 = l[i] + p1 = l[-1] + l[i] = f(p_1, p0, p1) + + return [ x+y for x,y in zip(gs,es) ] + + +def compute(data): + + d = [] + for e0, p0, a0, e1, p1, a1 in data: + x = (e1+p1)-(e0+p0) + w = 1./sqrt(p0**2 + p1**2) + bias = (a1-1.)*p1 - (a0-1.)*p0 + d.append( (x,w,bias) ) + + x = [] + target = (scipy.stats.norm.cdf(1.)-0.5)*2 # = 0.6827 + + print("| %2s | %8s | %8s | %8s | %8s | %8s |"%( "N", "DE", "+/-", "bias", "P(G)", "J")) + print("|----+----------+----------+----------+----------+----------|") + xmax = (0.,0.,0.,0.,0.,0,0.) + for i in range(len(data)-1): + jb, mu, sigma, p = jarque_bera( [ (x,w) for (x,w,bias) in d[i:] ] ) + bias = sum ( [ w * e for (_,w,e) in d[i:] ] ) / sum ( [ w for (_,w,_) in d[i:] ] ) + mu = (mu+0.5*bias) * to_eV + sigma = sigma * to_eV + bias = bias * to_eV + n = len(data[i:]) + beta = student(0.5*(1.+target/p) ,n) + err = sigma * beta + 0.5*abs(bias) + print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |"%( n, mu, err, bias, p, jb)) + if n < 3 : + continue + y = (err, p, mu, err, jb,n,bias) + if p > xmax[1]: xmax = y + if p < 0.8: + continue + x.append(y) + + x = sorted(x) + + print("|----+----------+----------+----------+----------+----------|") + if x != []: + xmax = x[0] + _, p, mu, err, jb, n, bias = xmax + beta = student(0.5*(1.+target/p),n) + print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |\n"%(n, mu, err, bias, p, jb)) + + return mu, err, bias, p + +filename = sys.argv[1] +print(filename) +if len(sys.argv) > 2: + state = int(sys.argv[2]) +else: + state = 1 +data = read_data(filename,state) +mu, err, bias, _ = compute(data) +print(" %s: %8.3f +/- %5.3f eV\n"%(filename, mu, err)) + +import numpy as np +A = np.array( [ [ data[-1][1], 1. ], + [ data[-2][1], 1. ] ] ) +B = np.array( [ [ data[-1][0] ], + [ data[-2][0] ] ] ) +E0 = np.linalg.solve(A,B)[1] +A = np.array( [ [ data[-1][4], 1. ], + [ data[-2][4], 1. ] ] ) +B = np.array( [ [ data[-1][3] ], + [ data[-2][3] ] ] ) +E1 = np.linalg.solve(A,B)[1] +average_2 = (E1-E0)*to_eV + +A = np.array( [ [ data[-1][1], 1. ], + [ data[-2][1], 1. ], + [ data[-3][1], 1. ] ] ) +B = np.array( [ [ data[-1][0] ], + [ data[-2][0] ], + [ data[-3][0] ] ] ) +E0 = np.linalg.lstsq(A,B,rcond=None)[0][1] +A = np.array( [ [ data[-1][4], 1. ], + [ data[-2][4], 1. ], + [ data[-3][4], 1. ] ] ) +B = np.array( [ [ data[-1][3] ], + [ data[-2][3] ], + [ data[-3][3] ] ] ) +E1 = np.linalg.lstsq(A,B,rcond=None)[0][1] +average_3 = (E1-E0)*to_eV + +exc = ((data[-1][3] + data[-1][4]) - (data[-1][0] + data[-1][1])) * to_eV +error_2 = abs(average_2 - average_3) +error_3 = abs(average_3 - exc) +print(" 2-3 points: %.3f +/- %.3f "% (average_3, error_2)) +print(" largest wf: %.3f +/- %.3f "%(average_3, error_3)) + + From 918839fbf6636c36ae9752b2d029ff5f31d157d3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 01:22:24 +0200 Subject: [PATCH 11/20] Added JSON in FCI_TC --- src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 10 +++++++++- src/fci_tc_bi/fci_tc_bi_ortho.irp.f | 6 ++++++ src/tc_bi_ortho/h_tc_s2_u0.irp.f | 20 ++++++++++---------- 3 files changed, 25 insertions(+), 11 deletions(-) diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 1f4fe849..284a1e2e 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -94,7 +94,15 @@ subroutine run_stochastic_cipsi call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection ! stop - N_iter += 1 + call print_summary(psi_energy_with_nucl_rep, & + pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) + + call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) + + call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) + call print_extrapolated_energy() + call print_mol_properties() + call write_cipsi_json(pt2_data,pt2_data_err) if (qp_stop()) exit diff --git a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f index 84ac8166..ed75c882 100644 --- a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f +++ b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f @@ -62,6 +62,7 @@ subroutine run_cipsi_tc endif endif ! --- + write(json_unit,json_array_open_fmt) 'fci_tc' if (do_pt2) then call run_stochastic_cipsi @@ -69,6 +70,11 @@ subroutine run_cipsi_tc call run_cipsi endif + write(json_unit,json_dict_uopen_fmt) + write(json_unit,json_dict_close_fmtx) + write(json_unit,json_array_close_fmtx) + call json_close + else PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks if(elec_alpha_num+elec_beta_num.ge.3)then diff --git a/src/tc_bi_ortho/h_tc_s2_u0.irp.f b/src/tc_bi_ortho/h_tc_s2_u0.irp.f index 30b0f273..b9b85a96 100644 --- a/src/tc_bi_ortho/h_tc_s2_u0.irp.f +++ b/src/tc_bi_ortho/h_tc_s2_u0.irp.f @@ -12,9 +12,9 @@ subroutine get_H_tc_s2_l0_r0(l_0,r_0,N_st,sze,energies, s2) ! istart, iend, ishift, istep are used in ZMQ parallelization. END_DOC integer, intent(in) :: N_st,sze - double precision, intent(in) :: l_0(sze,N_st), r_0(sze,N_st) + double precision, intent(inout) :: l_0(sze,N_st), r_0(sze,N_st) double precision, intent(out) :: energies(N_st), s2(N_st) - logical :: do_right + logical :: do_right integer :: istate double precision, allocatable :: s_0(:,:), v_0(:,:) double precision :: u_dot_v, norm @@ -40,7 +40,7 @@ subroutine H_tc_s2_u_0_opt(v_0,s_0,u_0,N_st,sze) END_DOC integer, intent(in) :: N_st,sze double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st) - logical :: do_right + logical :: do_right do_right = .True. call H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right) end @@ -57,7 +57,7 @@ subroutine H_tc_s2_dagger_u_0_opt(v_0,s_0,u_0,N_st,sze) END_DOC integer, intent(in) :: N_st,sze double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st) - logical :: do_right + logical :: do_right do_right = .False. call H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right) end @@ -77,7 +77,7 @@ subroutine H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right) END_DOC integer, intent(in) :: N_st,sze double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st) - logical, intent(in) :: do_right + logical, intent(in) :: do_right integer :: k double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t @@ -124,7 +124,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishi use bitmasks implicit none BEGIN_DOC - ! Computes $v_t = H | u_t\rangle$ + ! Computes $v_t = H | u_t\rangle$ ! ! Default should be 1,N_det,0,1 ! @@ -132,7 +132,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishi END_DOC integer, intent(in) :: N_st,sze,istart,iend,ishift,istep double precision, intent(in) :: u_t(N_st,N_det) - logical, intent(in) :: do_right + logical, intent(in) :: do_right double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze) @@ -165,7 +165,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,ie END_DOC integer, intent(in) :: N_st,sze,istart,iend,ishift,istep double precision, intent(in) :: u_t(N_st,N_det) - logical, intent(in) :: do_right + logical, intent(in) :: do_right double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze) double precision :: hij, sij @@ -542,7 +542,7 @@ compute_singles=.True. lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow) + tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow) ! call i_H_j( tmp_det, tmp_det2, $N_int, hij) ! call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij) if(do_right)then @@ -693,7 +693,7 @@ compute_singles=.True. lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol) + tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol) ! call i_H_j( tmp_det, tmp_det2, $N_int, hij) ! call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij) if(do_right)then From 54a88fe4caf88dc8be032d922fe87c0b250deab5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 01:32:05 +0200 Subject: [PATCH 12/20] Added JSON to fci_tc_bi --- src/cipsi_tc_bi_ortho/NEED | 3 +- src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 2 +- src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f | 53 ++++++++++ src/fci_tc_bi/NEED | 1 + src/fci_tc_bi/save_energy.irp.f | 4 +- src/iterations_tc/EZFIO.cfg | 24 ----- src/iterations_tc/NEED | 0 src/iterations_tc/io.irp.f | 37 ------- src/iterations_tc/iterations.irp.f | 43 -------- src/iterations_tc/print_extrapolation.irp.f | 46 -------- src/iterations_tc/print_summary.irp.f | 104 ------------------- 11 files changed, 59 insertions(+), 258 deletions(-) create mode 100644 src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f delete mode 100644 src/iterations_tc/EZFIO.cfg delete mode 100644 src/iterations_tc/NEED delete mode 100644 src/iterations_tc/io.irp.f delete mode 100644 src/iterations_tc/iterations.irp.f delete mode 100644 src/iterations_tc/print_extrapolation.irp.f delete mode 100644 src/iterations_tc/print_summary.irp.f diff --git a/src/cipsi_tc_bi_ortho/NEED b/src/cipsi_tc_bi_ortho/NEED index 4dd1af36..8f05be69 100644 --- a/src/cipsi_tc_bi_ortho/NEED +++ b/src/cipsi_tc_bi_ortho/NEED @@ -1,6 +1,7 @@ +json mpi perturbation zmq -iterations_tc +iterations csf tc_bi_ortho diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 284a1e2e..a06f28e9 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -101,7 +101,7 @@ subroutine run_stochastic_cipsi call increment_n_iter(psi_energy_with_nucl_rep, pt2_data) call print_extrapolated_energy() - call print_mol_properties() +! call print_mol_properties() call write_cipsi_json(pt2_data,pt2_data_err) if (qp_stop()) exit diff --git a/src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f b/src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f new file mode 100644 index 00000000..98a402a2 --- /dev/null +++ b/src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f @@ -0,0 +1,53 @@ +subroutine write_cipsi_json(pt2_data, pt2_data_err) + use selection_types + implicit none + BEGIN_DOC +! Writes JSON data for CIPSI runs + END_DOC + type(pt2_type), intent(in) :: pt2_data, pt2_data_err + integer :: i,j,k + + call lock_io + character*(64), allocatable :: fmtk(:) + integer :: N_states_p, N_iter_p + N_states_p = min(N_states,N_det) + N_iter_p = min(N_iter,8) + allocate(fmtk(0:N_iter_p)) + fmtk(:) = '('' '',E22.15,'','')' + fmtk(N_iter_p) = '('' '',E22.15)' + + write(json_unit, json_dict_uopen_fmt) + write(json_unit, json_int_fmt) 'n_det', N_det + if (s2_eig) then + write(json_unit, json_int_fmt) 'n_cfg', N_configuration + if (only_expected_s2) then + write(json_unit, json_int_fmt) 'n_csf', N_csf + endif + endif + write(json_unit, json_array_open_fmt) 'states' + do k=1,N_states_p + write(json_unit, json_dict_uopen_fmt) + write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k) + write(json_unit, json_real_fmt) 's2', psi_s2(k) + write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k) + write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k) + write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k) + write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k) + write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k) + write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k) + write(json_unit, json_array_open_fmt) 'ex_energy' + do i=2,N_iter_p + write(json_unit, fmtk(i)) extrapolated_energy(i,k) + enddo + write(json_unit, json_array_close_fmtx) + if (k < N_states_p) then + write(json_unit, json_dict_close_fmt) + else + write(json_unit, json_dict_close_fmtx) + endif + enddo + write(json_unit, json_array_close_fmtx) + write(json_unit, json_dict_close_fmt) + deallocate(fmtk) + call unlock_io +end diff --git a/src/fci_tc_bi/NEED b/src/fci_tc_bi/NEED index 000b0deb..3bb9515a 100644 --- a/src/fci_tc_bi/NEED +++ b/src/fci_tc_bi/NEED @@ -1,3 +1,4 @@ +json tc_bi_ortho davidson_undressed cipsi_tc_bi_ortho diff --git a/src/fci_tc_bi/save_energy.irp.f b/src/fci_tc_bi/save_energy.irp.f index 7c41d00f..421ae5f8 100644 --- a/src/fci_tc_bi/save_energy.irp.f +++ b/src/fci_tc_bi/save_energy.irp.f @@ -4,6 +4,6 @@ subroutine save_energy(E,pt2) ! Saves the energy in |EZFIO|. END_DOC double precision, intent(in) :: E(N_states), pt2(N_states) - call ezfio_set_fci_tc_energy(E(1:N_states)) - call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states)) + call ezfio_set_fci_tc_bi_energy(E(1:N_states)) + call ezfio_set_fci_tc_bi_energy_pt2(E(1:N_states)+pt2(1:N_states)) end diff --git a/src/iterations_tc/EZFIO.cfg b/src/iterations_tc/EZFIO.cfg deleted file mode 100644 index 2a5e94a7..00000000 --- a/src/iterations_tc/EZFIO.cfg +++ /dev/null @@ -1,24 +0,0 @@ -[n_iter] -interface: ezfio -doc: Number of saved iterations -type:integer -default: 1 - -[n_det_iterations] -interface: ezfio, provider -doc: Number of determinants at each iteration -type: integer -size: (100) - -[energy_iterations] -interface: ezfio, provider -doc: The variational energy at each iteration -type: double precision -size: (determinants.n_states,100) - -[pt2_iterations] -interface: ezfio, provider -doc: The |PT2| correction at each iteration -type: double precision -size: (determinants.n_states,100) - diff --git a/src/iterations_tc/NEED b/src/iterations_tc/NEED deleted file mode 100644 index e69de29b..00000000 diff --git a/src/iterations_tc/io.irp.f b/src/iterations_tc/io.irp.f deleted file mode 100644 index 821f5e84..00000000 --- a/src/iterations_tc/io.irp.f +++ /dev/null @@ -1,37 +0,0 @@ -BEGIN_PROVIDER [ integer, n_iter ] - implicit none - BEGIN_DOC -! number of iterations - END_DOC - - logical :: has - PROVIDE ezfio_filename - if (mpi_master) then - - double precision :: zeros(N_states,100) - integer :: izeros(100) - zeros = 0.d0 - izeros = 0 - call ezfio_set_iterations_n_iter(0) - call ezfio_set_iterations_energy_iterations(zeros) - call ezfio_set_iterations_pt2_iterations(zeros) - call ezfio_set_iterations_n_det_iterations(izeros) - n_iter = 1 - endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read n_iter with MPI' - endif - IRP_ENDIF - - call write_time(6) - -END_PROVIDER - diff --git a/src/iterations_tc/iterations.irp.f b/src/iterations_tc/iterations.irp.f deleted file mode 100644 index 2f1cf0c1..00000000 --- a/src/iterations_tc/iterations.irp.f +++ /dev/null @@ -1,43 +0,0 @@ -BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ] - implicit none - BEGIN_DOC - ! Extrapolated energy, using E_var = f(PT2) where PT2=0 - END_DOC -! integer :: i - extrapolated_energy = 0.D0 -END_PROVIDER - - subroutine get_extrapolated_energy(Niter,ept2,pt1,extrap_energy) - implicit none - integer, intent(in) :: Niter - double precision, intent(in) :: ept2(Niter),pt1(Niter),extrap_energy(Niter) - call extrapolate_data(Niter,ept2,pt1,extrap_energy) - end - -subroutine save_iterations(e_, pt2_,n_) - implicit none - BEGIN_DOC -! Update the energy in the EZFIO file. - END_DOC - integer, intent(in) :: n_ - double precision, intent(in) :: e_(N_states), pt2_(N_states) - integer :: i - - if (N_iter == 101) then - do i=2,N_iter-1 - energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter) - pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) - enddo - N_iter = N_iter-1 - TOUCH N_iter - endif - - energy_iterations(1:N_states,N_iter) = e_(1:N_states) - pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states) - n_det_iterations(N_iter) = n_ - call ezfio_set_iterations_N_iter(N_iter) - call ezfio_set_iterations_energy_iterations(energy_iterations) - call ezfio_set_iterations_pt2_iterations(pt2_iterations) - call ezfio_set_iterations_n_det_iterations(n_det_iterations) -end - diff --git a/src/iterations_tc/print_extrapolation.irp.f b/src/iterations_tc/print_extrapolation.irp.f deleted file mode 100644 index cb46fb67..00000000 --- a/src/iterations_tc/print_extrapolation.irp.f +++ /dev/null @@ -1,46 +0,0 @@ -subroutine print_extrapolated_energy - implicit none - BEGIN_DOC -! Print the extrapolated energy in the output - END_DOC - - integer :: i,k - - if (N_iter< 2) then - return - endif - write(*,'(A)') '' - write(*,'(A)') 'Extrapolated energies' - write(*,'(A)') '------------------------' - write(*,'(A)') '' - - print *, '' - print *, 'State ', 1 - print *, '' - write(*,*) '=========== ', '===================' - write(*,*) 'minimum PT2 ', 'Extrapolated energy' - write(*,*) '=========== ', '===================' - do k=2,min(N_iter,8) - write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1) - enddo - write(*,*) '=========== ', '===================' - - do i=2, min(N_states,N_det) - print *, '' - print *, 'State ', i - print *, '' - write(*,*) '=========== ', '=================== ', '=================== ', '===================' - write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) ' - write(*,*) '=========== ', '=================== ', '=================== ', '===================' - do k=2,min(N_iter,8) - write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), & - extrapolated_energy(k,i) - extrapolated_energy(k,1), & - (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0 - enddo - write(*,*) '=========== ', '=================== ', '=================== ', '===================' - enddo - - print *, '' - -end subroutine - diff --git a/src/iterations_tc/print_summary.irp.f b/src/iterations_tc/print_summary.irp.f deleted file mode 100644 index 8e6285e2..00000000 --- a/src/iterations_tc/print_summary.irp.f +++ /dev/null @@ -1,104 +0,0 @@ -subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_) - use selection_types - implicit none - BEGIN_DOC -! Print the extrapolated energy in the output - END_DOC - - integer, intent(in) :: n_det_, n_configuration_, n_st - double precision, intent(in) :: e_(n_st), s2_(n_st) - type(pt2_type) , intent(in) :: pt2_data, pt2_data_err - integer :: i, k - integer :: N_states_p - character*(9) :: pt2_string - character*(512) :: fmt - - if (do_pt2) then - pt2_string = ' ' - else - pt2_string = '(approx)' - endif - - N_states_p = min(N_det_,n_st) - - print *, '' - print '(A,I12)', 'Summary at N_det = ', N_det_ - print '(A)', '-----------------------------------' - print *, '' - - write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' - write(*,fmt) - write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))' - write(*,fmt) ('State',k, k=1,N_states_p) - write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' - write(*,fmt) - write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))' - write(*,fmt) '# E ', e_(1:N_states_p) - if (N_states_p > 1) then - write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) - write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 - endif - write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' - write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) - write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p) - write(*,'(A)') '#' - write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p) - write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p) - if (N_states_p > 1) then - write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & - dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) - write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & - dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p) - endif - write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' - write(*,fmt) - print *, '' - - print *, 'N_det = ', N_det_ - print *, 'N_states = ', n_st - if (s2_eig) then - print *, 'N_cfg = ', N_configuration_ - if (only_expected_s2) then - print *, 'N_csf = ', N_csf - endif - endif - print *, '' - - do k=1, N_states_p - print*,'* State ',k - print *, '< S^2 > = ', s2_(k) - print *, 'E = ', e_(k) - print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k) - print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k)) - print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) - print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) - print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) - print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) - print *, '' - enddo - - print *, '-----' - if(n_st.gt.1)then - print *, 'Variational Energy difference (au | eV)' - do i=2, N_states_p - print*,'Delta E = ', (e_(i) - e_(1)), & - (e_(i) - e_(1)) * 27.211396641308d0 - enddo - print *, '-----' - print*, 'Variational + perturbative Energy difference (au | eV)' - do i=2, N_states_p - print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), & - (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0 - enddo - print *, '-----' - print*, 'Variational + renormalized perturbative Energy difference (au | eV)' - do i=2, N_states_p - print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), & - (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0 - enddo - endif - -! call print_energy_components() - -end subroutine - From 67902d437776d5caeb4f9292796baafaa9ca66fe Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 01:36:49 +0200 Subject: [PATCH 13/20] Removed print_e_conv: should be replaced by python script --- src/tools/print_e_conv.irp.f | 80 ------------------------------------ 1 file changed, 80 deletions(-) delete mode 100644 src/tools/print_e_conv.irp.f diff --git a/src/tools/print_e_conv.irp.f b/src/tools/print_e_conv.irp.f deleted file mode 100644 index e358ebc1..00000000 --- a/src/tools/print_e_conv.irp.f +++ /dev/null @@ -1,80 +0,0 @@ -program print_e_conv - implicit none - BEGIN_DOC -! program that prints in a human readable format the convergence of the CIPSI algorithm. -! -! for all istate, this program produces -! -! * a file "EZFIO.istate.conv" containing the variational and var+PT2 energies as a function of N_det -! -! * for istate > 1, a file EZFIO.istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det - END_DOC - - provide ezfio_filename - call routine_e_conv - end - -subroutine routine_e_conv - implicit none - BEGIN_DOC -! routine called by :c:func:`print_e_conv` - END_DOC - integer :: N_iter_tmp - integer :: i,istate - character*(128) :: output - integer :: i_unit_output,getUnitAndOpen - character*(128) :: filename - - integer, allocatable :: n_det_tmp(:) - call ezfio_get_iterations_N_iter(N_iter_tmp) - print*,'N_iter_tmp = ',N_iter_tmp - double precision, allocatable :: e(:,:),pt2(:,:) - allocate(e(N_states, 100),pt2(N_states, 100),n_det_tmp(100)) - call ezfio_get_iterations_energy_iterations(e) - call ezfio_get_iterations_pt2_iterations(pt2) - call ezfio_get_iterations_n_det_iterations(n_det_tmp) - - - do istate = 1, N_states - if (istate.lt.10)then - write (filename, "(I1)")istate - else - write (filename, "(I2)")istate - endif - print*,filename - output=trim(ezfio_filename)//'.'//trim(filename)//'.conv' - output=trim(output) - print*,'output = ',trim(output) - i_unit_output = getUnitAndOpen(output,'w') - write(i_unit_output,*)'# N_det E_var E_var + PT2' - do i = 1, N_iter_tmp - write(i_unit_output,'(I9,X,3(F16.10,X))')n_det_tmp(i),e(istate,i),e(istate,i) + pt2(istate,i) - enddo - enddo - - if(N_states.gt.1)then - double precision, allocatable :: deltae(:,:),deltae_pt2(:,:) - allocate(deltae(N_states,100),deltae_pt2(N_states,100)) - do i = 1, N_iter_tmp - do istate = 1, N_states - deltae(istate,i) = e(istate,i) - e(1,i) - deltae_pt2(istate,i) = e(istate,i) + pt2(istate,i) - (e(1,i) + pt2(1,i)) - enddo - enddo - do istate = 2, N_states - if (istate.lt.10)then - write (filename, "(I1)")istate - else - write (filename, "(I2)")istate - endif - output=trim(ezfio_filename)//'.'//trim(filename)//'.delta_e.conv' - print*,'output = ',trim(output) - i_unit_output = getUnitAndOpen(output,'w') - write(i_unit_output,*)'# N_det Delta E_var Delta (E_var + PT2)' - do i = 1, N_iter_tmp - write(i_unit_output,'(I9,X,100(F16.10,X))')n_det_tmp(i),deltae(istate,i),deltae_pt2(istate,i) - enddo - enddo - endif - -end From 64bfddbb00e50d46bbae11c2cc895435d298bcd0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 10:48:21 +0200 Subject: [PATCH 14/20] Renamed scripts/exc_energy_error.py scripts/qp_exc_energy.py --- scripts/{exc_energy_error.py => qp_exc_energy.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/{exc_energy_error.py => qp_exc_energy.py} (100%) diff --git a/scripts/exc_energy_error.py b/scripts/qp_exc_energy.py similarity index 100% rename from scripts/exc_energy_error.py rename to scripts/qp_exc_energy.py From d3b76e5957aadc11f2090084491e8b90bb6a50ff Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 14 Apr 2023 10:56:07 +0200 Subject: [PATCH 15/20] changed h_p to h --- external/qp2-dependencies | 2 +- src/davidson/diagonalization_hs2_dressed.irp.f | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/external/qp2-dependencies b/external/qp2-dependencies index fd43778e..e0d0e02e 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit fd43778e12bb5858c4c780c34346be0f158b8cc7 +Subproject commit e0d0e02e9f5ece138d1520106954a881ab0b8db2 diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 8117f320..d37b7386 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -465,8 +465,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ integer :: lwork, info double precision, allocatable :: work(:) -! y = h - y = h_p + y = h lwork = -1 allocate(work(1)) call dsygv(1,'V','U',shift2,y,size(y,1), & From c80ebe27b8d2d8592570aecb46e19fac6ca65064 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Apr 2023 10:31:24 +0200 Subject: [PATCH 16/20] Introducing Cholesky-decomposed SCF --- src/ao_two_e_ints/cholesky.irp.f | 100 +++++++++ src/ao_two_e_ints/map_integrals.irp.f | 2 +- src/hartree_fock/fock_matrix_hf.irp.f | 311 ++++++++++++++++---------- src/utils/linear_algebra.irp.f | 2 +- 4 files changed, 298 insertions(+), 117 deletions(-) create mode 100644 src/ao_two_e_ints/cholesky.irp.f diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f new file mode 100644 index 00000000..d4c201aa --- /dev/null +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -0,0 +1,100 @@ +BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ] + implicit none + BEGIN_DOC + ! Number of Cholesky vectors in AO basis + END_DOC + + integer :: i,j,k,l + double precision :: xnorm0, x, integral + double precision, external :: ao_two_e_integral + + cholesky_ao_num_guess = 0 + xnorm0 = 0.d0 + x = 0.d0 + do j=1,ao_num + do i=1,ao_num + integral = ao_two_e_integral(i,i,j,j) + if (integral > ao_integrals_threshold) then + cholesky_ao_num_guess += 1 + else + x += integral + endif + enddo + enddo + print *, 'Cholesky decomposition of AO integrals' + print *, '--------------------------------------' + print *, '' + print *, 'Estimated Error: ', x + print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)' + +END_PROVIDER + + BEGIN_PROVIDER [ integer, cholesky_ao_num ] +&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, cholesky_ao_num_guess) ] + use mmap_module + implicit none + BEGIN_DOC + ! Cholesky vectors in AO basis: (ik|a): + ! = (ik|jl) = sum_a (ik|a).(a|jl) + END_DOC + + type(c_ptr) :: ptr + integer :: fd, i,j,k,l, rank + double precision, pointer :: ao_integrals(:,:,:,:) + double precision, external :: ao_two_e_integral + + ! Store AO integrals in a memory mapped file + call mmap(trim(ezfio_work_dir)//'ao_integrals', & + (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), & + 8, fd, .False., ptr) + call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/)) + + double precision :: integral + logical, external :: ao_two_e_integral_zero + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic) + do l=1,ao_num + do j=1,l + do k=1,ao_num + do i=1,k + if (ao_two_e_integral_zero(i,j,k,l)) cycle + integral = ao_two_e_integral(i,k,j,l) + 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 + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + ! Call Lapack + cholesky_ao_num = cholesky_ao_num_guess + call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao) + print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)' + + ! Remove mmap + double precision, external :: getUnitAndOpen + call munmap( & + (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), & + 8, fd, ptr) + open(unit=99,file=trim(ezfio_work_dir)//'ao_integrals') + close(99, status='delete') + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ] + implicit none + BEGIN_DOC +! Transposed of the Cholesky vectors in AO basis set + END_DOC + integer :: i,j,k + do j=1,ao_num + do i=1,ao_num + do k=1,ao_num + cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index fa7c29cc..7d6a7da4 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -486,7 +486,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) PROVIDE ao_two_e_integrals_in_map ao_integrals_map if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 + out_val(1:sze) = 0.d0 return endif diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index d7d8fa7d..12641516 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -15,115 +15,59 @@ double precision, allocatable :: ao_two_e_integral_alpha_tmp(:,:) double precision, allocatable :: ao_two_e_integral_beta_tmp(:,:) - ao_two_e_integral_alpha = 0.d0 - ao_two_e_integral_beta = 0.d0 - if (do_direct_integrals) then + if (.True.) then ! Use Cholesky-decomposed integrals + ao_two_e_integral_alpha(:,:) = ao_two_e_integral_alpha_chol(:,:) + ao_two_e_integral_beta (:,:) = ao_two_e_integral_beta_chol (:,:) - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0, & - !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2, & - !$OMP local_threshold)& - !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& - !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, & - !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta) + else ! Use integrals in AO basis set - allocate(keys(1), values(1)) - allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & - ao_two_e_integral_beta_tmp(ao_num,ao_num)) - ao_two_e_integral_alpha_tmp = 0.d0 - ao_two_e_integral_beta_tmp = 0.d0 + ao_two_e_integral_alpha = 0.d0 + ao_two_e_integral_beta = 0.d0 + if (do_direct_integrals) then - q = ao_num*ao_num*ao_num*ao_num - !$OMP DO SCHEDULE(static,64) - do p=1_8,q - call two_e_integrals_index_reverse(kk,ii,ll,jj,p) - if ( (kk(1)>ao_num).or. & - (ii(1)>ao_num).or. & - (jj(1)>ao_num).or. & - (ll(1)>ao_num) ) then - cycle - endif - k = kk(1) - i = ii(1) - l = ll(1) - j = jj(1) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,& + !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,& + !$OMP local_threshold) & + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,& + !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta) - logical, external :: ao_two_e_integral_zero - if (ao_two_e_integral_zero(i,k,j,l)) then - cycle - endif - local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) - if (local_threshold < ao_integrals_threshold) then - cycle - endif - i0 = i - j0 = j - k0 = k - l0 = l - values(1) = 0.d0 - local_threshold = ao_integrals_threshold/local_threshold - do k2=1,8 - if (kk(k2)==0) then - cycle - endif - i = ii(k2) - j = jj(k2) - k = kk(k2) - l = ll(k2) - c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) - c1 = SCF_density_matrix_ao_alpha(k,i) - c2 = SCF_density_matrix_ao_beta(k,i) - if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then - cycle - endif - if (values(1) == 0.d0) then - values(1) = ao_two_e_integral(k0,l0,i0,j0) - endif - integral = c0 * values(1) - ao_two_e_integral_alpha_tmp(i,j) += integral - ao_two_e_integral_beta_tmp (i,j) += integral - integral = values(1) - ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral - ao_two_e_integral_beta_tmp (l,j) -= c2 * integral - enddo - enddo - !$OMP END DO NOWAIT - !$OMP CRITICAL - ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp - ao_two_e_integral_beta += ao_two_e_integral_beta_tmp - !$OMP END CRITICAL - deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) - !$OMP END PARALLEL - else - PROVIDE ao_two_e_integrals_in_map + allocate(keys(1), values(1)) + allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_num)) + ao_two_e_integral_alpha_tmp = 0.d0 + ao_two_e_integral_beta_tmp = 0.d0 - integer(omp_lock_kind) :: lck(ao_num) - integer(map_size_kind) :: i8 - integer :: ii(8), jj(8), kk(8), ll(8), k2 - integer(cache_map_size_kind) :: n_elements_max, n_elements - integer(key_kind), allocatable :: keys(:) - double precision, allocatable :: values(:) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & - !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)& - !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& - !$OMP ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta) - - call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) - allocate(keys(n_elements_max), values(n_elements_max)) - allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & - ao_two_e_integral_beta_tmp(ao_num,ao_num)) - ao_two_e_integral_alpha_tmp = 0.d0 - ao_two_e_integral_beta_tmp = 0.d0 - - !$OMP DO SCHEDULE(static,1) - do i8=0_8,ao_integrals_map%map_size - n_elements = n_elements_max - call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) - do k1=1,n_elements - call two_e_integrals_index_reverse(kk,ii,ll,jj,keys(k1)) + q = ao_num*ao_num*ao_num*ao_num + !$OMP DO SCHEDULE(static,64) + do p=1_8,q + call two_e_integrals_index_reverse(kk,ii,ll,jj,p) + if ( (kk(1)>ao_num).or. & + (ii(1)>ao_num).or. & + (jj(1)>ao_num).or. & + (ll(1)>ao_num) ) then + cycle + endif + k = kk(1) + i = ii(1) + l = ll(1) + j = jj(1) + logical, external :: ao_two_e_integral_zero + if (ao_two_e_integral_zero(i,k,j,l)) then + cycle + endif + local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) + if (local_threshold < ao_integrals_threshold) then + cycle + endif + i0 = i + j0 = j + k0 = k + l0 = l + values(1) = 0.d0 + local_threshold = ao_integrals_threshold/local_threshold do k2=1,8 if (kk(k2)==0) then cycle @@ -132,25 +76,162 @@ j = jj(k2) k = kk(k2) l = ll(k2) - integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1) + c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) + c1 = SCF_density_matrix_ao_alpha(k,i) + c2 = SCF_density_matrix_ao_beta(k,i) + if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then + cycle + endif + if (values(1) == 0.d0) then + values(1) = ao_two_e_integral(k0,l0,i0,j0) + endif + integral = c0 * values(1) ao_two_e_integral_alpha_tmp(i,j) += integral ao_two_e_integral_beta_tmp (i,j) += integral - integral = values(k1) - ao_two_e_integral_alpha_tmp(l,j) -= SCF_density_matrix_ao_alpha(k,i) * integral - ao_two_e_integral_beta_tmp (l,j) -= SCF_density_matrix_ao_beta (k,i) * integral + integral = values(1) + ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral + ao_two_e_integral_beta_tmp (l,j) -= c2 * integral enddo enddo - enddo - !$OMP END DO NOWAIT - !$OMP CRITICAL - ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp - ao_two_e_integral_beta += ao_two_e_integral_beta_tmp - !$OMP END CRITICAL - deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) - !$OMP END PARALLEL + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + else + PROVIDE ao_two_e_integrals_in_map + integer(omp_lock_kind) :: lck(ao_num) + integer(map_size_kind) :: i8 + integer :: ii(8), jj(8), kk(8), ll(8), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max,& + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_num)) + ao_two_e_integral_alpha_tmp = 0.d0 + ao_two_e_integral_beta_tmp = 0.d0 + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + call two_e_integrals_index_reverse(kk,ii,ll,jj,keys(k1)) + + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1) + ao_two_e_integral_alpha_tmp(i,j) += integral + ao_two_e_integral_beta_tmp (i,j) += integral + integral = values(k1) + ao_two_e_integral_alpha_tmp(l,j) -= SCF_density_matrix_ao_alpha(k,i) * integral + ao_two_e_integral_beta_tmp (l,j) -= SCF_density_matrix_ao_beta (k,i) * integral + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + endif endif +END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha_chol, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta_chol , (ao_num, ao_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha and Beta Fock matrices in AO basis set + END_DOC + + integer :: m,n,l,s,j + double precision :: integral + double precision, allocatable :: X(:), X2(:,:,:,:), X3(:,:,:,:) + + allocate (X(cholesky_ao_num)) + + + ! X(j) = \sum_{mn} SCF_density_matrix_ao(m,n) * cholesky_ao(m,n,j) + call dgemm('T','N',cholesky_ao_num,1,ao_num*ao_num,1.d0, & + cholesky_ao, ao_num*ao_num, & + SCF_density_matrix_ao, ao_num*ao_num,0.d0, & + X, cholesky_ao_num) +! + + ! ao_two_e_integral_alpha(m,n) = \sum_{j} cholesky_ao(m,n,j) * X(j) + call dgemm('N','N',ao_num*ao_num,1,cholesky_ao_num, 1.d0, & + cholesky_ao, ao_num*ao_num, & + X, cholesky_ao_num, 0.d0, & + ao_two_e_integral_alpha_chol, ao_num*ao_num) + + deallocate(X) + + ao_two_e_integral_beta_chol = ao_two_e_integral_alpha_chol + + + allocate(X2(ao_num,ao_num,cholesky_ao_num,2)) + +! ao_two_e_integral_alpha_chol (l,s) -= cholesky_ao(l,m,j) * SCF_density_matrix_ao_beta (m,n) * cholesky_ao(n,s,j) + + call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0, & + SCF_density_matrix_ao_alpha, ao_num, & + cholesky_ao, ao_num, 0.d0, & + X2(1,1,1,1), ao_num) + + call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0, & + SCF_density_matrix_ao_beta, ao_num, & + cholesky_ao, ao_num, 0.d0, & + X2(1,1,1,2), ao_num) + + allocate(X3(ao_num,cholesky_ao_num,ao_num,2)) + + do s=1,ao_num + do j=1,cholesky_ao_num + do m=1,ao_num + X3(m,j,s,1) = X2(m,s,j,1) + X3(m,j,s,2) = X2(m,s,j,2) + enddo + enddo + enddo + + deallocate(X2) + + call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0, & + cholesky_ao, ao_num, & + X3(1,1,1,1), ao_num*cholesky_ao_num, 1.d0, & + ao_two_e_integral_alpha_chol, ao_num) + + call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0, & + cholesky_ao, ao_num, & + X3(1,1,1,2), ao_num*cholesky_ao_num, 1.d0, & + ao_two_e_integral_beta_chol, ao_num) + + deallocate(X3) + END_PROVIDER BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index c02560e3..3b43d607 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1854,7 +1854,7 @@ do k = 1, N 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 +U(:,:) = 0.00D+0 do k = 1, N l = piv(k) U(l, :) = A(1:rank, k) From bfdfb546bd29d6d16d3cb2024f707f758c89a68f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 May 2023 11:41:17 +0200 Subject: [PATCH 17/20] Fix pt2_max extra iterations --- src/cipsi/cipsi.irp.f | 8 +------- src/cipsi/stochastic_cipsi.irp.f | 8 +------- 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 88aaeae0..1279b5cd 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -128,13 +128,7 @@ subroutine run_cipsi if (qp_stop()) exit enddo - if (.not.qp_stop()) then - if (N_det < N_det_max) then - call diagonalize_CI - call save_wavefunction - call save_energy(psi_energy_with_nucl_rep, zeros) - endif - + if ((.not.qp_stop()).and.(N_det > N_det_max)) then if (do_pt2) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index b83e658a..deffaeb7 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -118,13 +118,7 @@ subroutine run_stochastic_cipsi if (qp_stop()) exit enddo - if (.not.qp_stop()) then - if (N_det < N_det_max) then - call diagonalize_CI - call save_wavefunction - call save_energy(psi_energy_with_nucl_rep, zeros) - endif - + if ((.not.qp_stop()).and.(N_det > N_det_max)) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) From ba2e783e8c6531ddeda40ac72d223a1d1ee67dcc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 May 2023 11:41:17 +0200 Subject: [PATCH 18/20] Fix pt2_max extra iterations --- src/cipsi/cipsi.irp.f | 8 +------- src/cipsi/stochastic_cipsi.irp.f | 8 +------- 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 5225c6df..3b344f62 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -130,13 +130,7 @@ subroutine run_cipsi if (qp_stop()) exit enddo - if (.not.qp_stop()) then - if (N_det < N_det_max) then - call diagonalize_CI - call save_wavefunction - call save_energy(psi_energy_with_nucl_rep, zeros) - endif - + if ((.not.qp_stop()).and.(N_det > N_det_max)) then if (do_pt2) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 35e80eb8..6a14fbf3 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -119,13 +119,7 @@ subroutine run_stochastic_cipsi if (qp_stop()) exit enddo - if (.not.qp_stop()) then - if (N_det < N_det_max) then - call diagonalize_CI - call save_wavefunction - call save_energy(psi_energy_with_nucl_rep, zeros) - endif - + if ((.not.qp_stop()).and.(N_det > N_det_max)) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) From e8782546c1387bfc030f5e659d6a914d038d280f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 May 2023 11:54:29 +0200 Subject: [PATCH 19/20] Better fix for pt2_max --- src/cipsi/cipsi.irp.f | 8 +++++++- src/cipsi/stochastic_cipsi.irp.f | 8 +++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 3b344f62..f3a77609 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -130,7 +130,13 @@ subroutine run_cipsi if (qp_stop()) exit enddo - if ((.not.qp_stop()).and.(N_det > N_det_max)) then + ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2 + if ((.not.qp_stop()).and. & + (N_det > N_det_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.& + (correlation_energy_ratio <= correlation_energy_ratio_max) & + ) then if (do_pt2) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 6a14fbf3..0b9a3548 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -119,7 +119,13 @@ subroutine run_stochastic_cipsi if (qp_stop()) exit enddo - if ((.not.qp_stop()).and.(N_det > N_det_max)) then + ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2 + if ((.not.qp_stop()).and. & + (N_det > N_det_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.& + (correlation_energy_ratio <= correlation_energy_ratio_max) & + ) then call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) From 0330ac6ccc4eb24deaef70632aa9dcd1e8fabfaa Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 May 2023 15:50:40 +0200 Subject: [PATCH 20/20] 4idx transformation with cholesky --- config/ifort_2021_xHost.cfg | 2 +- src/ao_two_e_ints/EZFIO.cfg | 5 ++ src/hartree_fock/fock_matrix_hf.irp.f | 2 +- src/mo_two_e_ints/cholesky.irp.f | 16 ++++++ src/mo_two_e_ints/mo_bi_integrals.irp.f | 74 +++++++++++++++++++++++-- 5 files changed, 91 insertions(+), 8 deletions(-) create mode 100644 src/mo_two_e_ints/cholesky.irp.f diff --git a/config/ifort_2021_xHost.cfg b/config/ifort_2021_xHost.cfg index 1161833b..9170b059 100644 --- a/config/ifort_2021_xHost.cfg +++ b/config/ifort_2021_xHost.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable 5462 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=64 -DINTEL diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index b18c65d1..caed4698 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -18,3 +18,8 @@ interface: ezfio,provider,ocaml default: False ezfio_name: direct +[do_ao_cholesky] +type: logical +doc: Perform Cholesky decomposition of AO integrals +interface: ezfio,provider,ocaml +default: True diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index 12641516..8c6658c5 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -15,7 +15,7 @@ double precision, allocatable :: ao_two_e_integral_alpha_tmp(:,:) double precision, allocatable :: ao_two_e_integral_beta_tmp(:,:) - if (.True.) then ! Use Cholesky-decomposed integrals + if (do_ao_cholesky) then ! Use Cholesky-decomposed integrals ao_two_e_integral_alpha(:,:) = ao_two_e_integral_alpha_chol(:,:) ao_two_e_integral_beta (:,:) = ao_two_e_integral_beta_chol (:,:) diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f new file mode 100644 index 00000000..14d3c696 --- /dev/null +++ b/src/mo_two_e_ints/cholesky.irp.f @@ -0,0 +1,16 @@ +BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_ao_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in MO basis + END_DOC + + integer :: k + + !$OMP PARALLEL DO PRIVATE(k) + do k=1,cholesky_ao_num + call ao_to_mo(cholesky_ao(1,1,k),ao_num,cholesky_mo(1,1,k),mo_num) + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index ae299e9f..b7ef901d 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -50,13 +50,16 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] call cpu_time(cpu_1) if(no_vvvv_integrals)then -! call four_idx_novvvv call four_idx_novvvv_old else - if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then - call four_idx_dgemm + if (do_ao_cholesky) then + call add_integrals_to_map_cholesky else - call add_integrals_to_map(full_ijkl_bitmask_4) + if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then + call four_idx_dgemm + else + call add_integrals_to_map(full_ijkl_bitmask_4) + endif endif endif @@ -175,7 +178,7 @@ subroutine add_integrals_to_map(mask_ijkl) implicit none BEGIN_DOC - ! Adds integrals to tha MO map according to some bitmask + ! Adds integrals to the MO map according to some bitmask END_DOC integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) @@ -450,13 +453,72 @@ subroutine add_integrals_to_map(mask_ijkl) end +subroutine add_integrals_to_map_cholesky + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to the MO map using Cholesky vectors + END_DOC + + integer :: i,j,k,l,m + integer :: size_buffer, n_integrals + size_buffer = min(mo_num*mo_num*mo_num,16000000) + + double precision, allocatable :: Vtmp(:,:,:,:) + integer(key_kind) , allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_value(:) + + if (.True.) then + ! In-memory transformation + + allocate (Vtmp(mo_num,mo_num,mo_num,mo_num)) + + call dgemm('N','T',mo_num*mo_num,mo_num*mo_num,cholesky_ao_num,1.d0, & + cholesky_mo, mo_num*mo_num, & + cholesky_mo, mo_num*mo_num, 0.d0, & + Vtmp, mo_num*mo_num) + + !$OMP PARALLEL PRIVATE(i,j,k,l,n_integrals,buffer_value, buffer_i) + allocate (buffer_i(size_buffer), buffer_value(size_buffer)) + n_integrals = 0 + !$OMP DO + do l=1,mo_num + do k=1,l + do j=1,mo_num + do i=1,j + if (abs(Vtmp(i,j,k,l)) > mo_integrals_threshold) then + n_integrals += 1 + buffer_value(n_integrals) = Vtmp(i,j,k,l) + !DIR$ FORCEINLINE + call mo_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + n_integrals = 0 + endif + endif + enddo + enddo + enddo + enddo + !$OMP END DO + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + + deallocate(Vtmp) + call map_unique(mo_integrals_map) + + endif + +end subroutine add_integrals_to_map_three_indices(mask_ijk) use bitmasks implicit none BEGIN_DOC - ! Adds integrals to tha MO map according to some bitmask + ! Adds integrals to the MO map according to some bitmask END_DOC integer(bit_kind), intent(in) :: mask_ijk(N_int,3)