From 4c0de615fb0a572212c2c3940c601f64e8cd0164 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Apr 2023 19:36:39 +0200 Subject: [PATCH 01/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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 4b1d384fb92a14437d0eebc42bd2bde651024111 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 2 May 2023 19:01:25 +0200 Subject: [PATCH 17/28] added Slater-type envelope --- src/non_h_ints_mu/jast_deriv.irp.f | 41 +++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index 31856a3d..3a06196c 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -157,7 +157,18 @@ double precision function j1b_nucl(r) integer :: i double precision :: a, d, e, x, y, z - if(j1b_type .eq. 103) then + if(j1b_type .eq. 102) then + + j1b_nucl = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d)) + enddo + + elseif(j1b_type .eq. 103) then j1b_nucl = 1.d0 do i = 1, nucl_num @@ -215,7 +226,29 @@ subroutine grad1_j1b_nucl(r, grad) double precision :: fact_x, fact_y, fact_z double precision :: ax_der, ay_der, az_der, a_expo - if(j1b_type .eq. 103) then + if(j1b_type .eq. 102) then + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = dsqrt(x*x + y*y + z*z) + e = a * dexp(-a*d) / d + + fact_x += e * x + fact_y += e * y + fact_z += e * z + enddo + + grad(1) = fact_x + grad(2) = fact_y + grad(3) = fact_z + + elseif(j1b_type .eq. 103) then x = r(1) y = r(2) @@ -254,7 +287,7 @@ subroutine grad1_j1b_nucl(r, grad) grad(2) = fact_y grad(3) = fact_z - else if(j1b_type .eq. 104) then + elseif(j1b_type .eq. 104) then fact_x = 0.d0 fact_y = 0.d0 @@ -276,7 +309,7 @@ subroutine grad1_j1b_nucl(r, grad) grad(2) = 2.d0 * fact_y grad(3) = 2.d0 * fact_z - else if(j1b_type .eq. 105) then + elseif(j1b_type .eq. 105) then fact_x = 0.d0 fact_y = 0.d0 From f4aec2957d06f2d155adcd979f77b9566866d0a8 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Tue, 2 May 2023 21:47:29 +0200 Subject: [PATCH 18/28] norm modif in print_tc_wf --- src/tc_bi_ortho/print_tc_wf.irp.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tc_bi_ortho/print_tc_wf.irp.f b/src/tc_bi_ortho/print_tc_wf.irp.f index 58a733a7..0cf3ca87 100644 --- a/src/tc_bi_ortho/print_tc_wf.irp.f +++ b/src/tc_bi_ortho/print_tc_wf.irp.f @@ -26,7 +26,8 @@ subroutine write_l_r_wf integer :: i print*,'Writing the left-right wf' do i = 1, N_det - write(i_unit_output,*)i,psi_l_coef_sorted_bi_ortho_left(i),psi_r_coef_sorted_bi_ortho_right(i) + write(i_unit_output,*)i, psi_l_coef_sorted_bi_ortho_left(i)/psi_l_coef_sorted_bi_ortho_left(1) & + , psi_r_coef_sorted_bi_ortho_right(i)/psi_r_coef_sorted_bi_ortho_right(1) enddo From 5eb3b69873f0b7bdacfcea72cbe8364a51c46bfb Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Thu, 4 May 2023 01:42:06 +0200 Subject: [PATCH 19/28] mu(r) added --- src/non_h_ints_mu/jast_deriv.irp.f | 103 +++++++++++++++++++++++++++-- src/non_h_ints_mu/tc_integ.irp.f | 62 ++++++++--------- src/tc_keywords/EZFIO.cfg | 6 ++ 3 files changed, 134 insertions(+), 37 deletions(-) diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index 3a06196c..0adc49fb 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -74,6 +74,42 @@ !$OMP END DO !$OMP END PARALLEL + elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, jpoint, r1, r2, grad1_u2b, dx, dy, dz) & + !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, & + !$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + call grad1_j12_mu(r1, r2, grad1_u2b) + + dx = grad1_u2b(1) + dy = grad1_u2b(2) + dz = grad1_u2b(3) + + grad1_u12_num(jpoint,ipoint,1) = dx + grad1_u12_num(jpoint,ipoint,2) = dy + grad1_u12_num(jpoint,ipoint,3) = dz + + grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + else print *, ' j1b_type = ', j1b_type, 'not implemented yet' @@ -91,16 +127,16 @@ double precision function j12_mu(r1, r2) implicit none double precision, intent(in) :: r1(3), r2(3) - double precision :: mu_r12, r12 + double precision :: mu_tmp, r12 if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) - mu_r12 = mu_erf * r12 + mu_tmp = mu_erf * r12 - j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf + j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf else @@ -116,6 +152,8 @@ end function j12_mu subroutine grad1_j12_mu(r1, r2, grad) + include 'constants.include.F' + implicit none double precision, intent(in) :: r1(3), r2(3) double precision, intent(out) :: grad(3) @@ -129,7 +167,7 @@ subroutine grad1_j12_mu(r1, r2, grad) dy = r1(2) - r2(2) dz = r1(3) - r2(3) - r12 = dsqrt( dx * dx + dy * dy + dz * dz ) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) if(r12 .lt. 1d-10) return tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 @@ -138,6 +176,28 @@ subroutine grad1_j12_mu(r1, r2, grad) grad(2) = tmp * dy grad(3) = tmp * dz + elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then + + double precision :: mu_val, mu_tmp, mu_der(3) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + + call mu_r_val_and_grad(r1, r2, mu_val, mu_der) + mu_tmp = mu_val * r12 + tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) + grad(1) = tmp * mu_der(1) + grad(2) = tmp * mu_der(2) + grad(3) = tmp * mu_der(3) + + if(r12 .lt. 1d-10) return + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 + grad(1) = grad(1) + tmp * dx + grad(2) = grad(2) + tmp * dy + grad(3) = grad(3) + tmp * dz + else print *, ' j1b_type = ', j1b_type, 'not implemented yet' @@ -343,3 +403,38 @@ end subroutine grad1_j1b_nucl ! --- +subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: mu_val, mu_der(3) + + if(j1b_type .eq. 200) then + + double precision :: r(3), dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) + double precision :: dm_tot, tmp + + PROVIDE mu_r_ct + r(1) = 0.5d0 * (r1(1) + r2(1)) + r(2) = 0.5d0 * (r1(2) + r2(2)) + r(3) = 0.5d0 * (r1(3) + r2(3)) + + call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) + dm_tot = dm_a(1) + dm_b(1) + mu_val = mu_r_ct * dsqrt(dm_tot) + tmp = 0.25d0 * mu_r_ct / dm_tot + mu_der(1) = tmp * (grad_dm_a(1,1) + grad_dm_b(1,1)) + mu_der(2) = tmp * (grad_dm_a(2,1) + grad_dm_b(2,1)) + mu_der(3) = tmp * (grad_dm_a(3,1) + grad_dm_b(3,1)) + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + return +end subroutine mu_r_val_and_grad + +! --- diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index f725d134..784947b4 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -35,13 +35,13 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f PROVIDE j1b_type if(read_tc_integ) then + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read") read(11) int2_grad1_u12_ao - endif - if(j1b_type .eq. 3) then + else - if(.not.read_tc_integ) then + if(j1b_type .eq. 3) then PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b @@ -73,32 +73,29 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f !$OMP END DO !$OMP END PARALLEL - endif + elseif(j1b_type .ge. 100) then - elseif(j1b_type .ge. 100) then + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + PROVIDE grad1_u12_num - PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - PROVIDE grad1_u12_num - - double precision, allocatable :: tmp(:,:,:) - allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) - tmp = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + double precision, allocatable :: tmp(:,:,:) + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) + tmp = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO SCHEDULE (static) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL - if(.not.read_tc_integ) then int2_grad1_u12_ao = 0.d0 do m = 1, 3 !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 & @@ -108,7 +105,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) enddo - !! these dgemm are equivalen to + !! these dgemm are equivalent to !!$OMP PARALLEL & !!$OMP DEFAULT (NONE) & !!$OMP PRIVATE (j, i, ipoint, jpoint, w) & @@ -132,15 +129,14 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f !enddo !!$OMP END DO !!$OMP END PARALLEL + + deallocate(tmp) + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + endif - - deallocate(tmp) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - endif if(write_tc_integ.and.mpi_master) then diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 85c8dac3..2e14488e 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -124,6 +124,12 @@ doc: type of 1-body Jastrow interface: ezfio, provider, ocaml default: 0 +[mu_r_ct] +type: double precision +doc: a parameter used to define mu(r) +interface: ezfio, provider, ocaml +default: 6.203504908994001e-1 + [thr_degen_tc] type: Threshold doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue From bfdfb546bd29d6d16d3cb2024f707f758c89a68f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 May 2023 11:41:17 +0200 Subject: [PATCH 20/28] 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 21/28] 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 22/28] 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 23/28] 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) From 868d5a162b6101bd897b5f659a06009d3115b17f Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Sat, 6 May 2023 18:25:06 +0200 Subject: [PATCH 24/28] jast 201 added --- src/mo_basis/EZFIO.cfg | 6 +++ src/mo_basis/mos_aux.irp.f | 53 ++++++++++++++++++++++ src/non_h_ints_mu/jast_deriv.irp.f | 65 +++++++++++++++++++++++---- src/tc_bi_ortho/print_tc_energy.irp.f | 4 ++ src/tc_keywords/EZFIO.cfg | 2 +- 5 files changed, 120 insertions(+), 10 deletions(-) create mode 100644 src/mo_basis/mos_aux.irp.f diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index 81ffba5c..4c4f1eca 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO| interface: ezfio size: (ao_basis.ao_num,mo_basis.mo_num) +[mo_coef_aux] +type: double precision +doc: AUX Coefficient of the i-th |AO| on the j-th |MO| +interface: ezfio +size: (ao_basis.ao_num,mo_basis.mo_num) + [mo_coef_imag] type: double precision doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO| diff --git a/src/mo_basis/mos_aux.irp.f b/src/mo_basis/mos_aux.irp.f new file mode 100644 index 00000000..27a874b1 --- /dev/null +++ b/src/mo_basis/mos_aux.irp.f @@ -0,0 +1,53 @@ + +! --- + +BEGIN_PROVIDER [double precision, mo_coef_aux, (ao_num,mo_num)] + + implicit none + integer :: i, j + logical :: exists + double precision, allocatable :: buffer(:,:) + + PROVIDE ezfio_filename + + if (mpi_master) then + ! Coefs + call ezfio_has_mo_basis_mo_coef_aux(exists) + 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(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_aux with MPI' + endif + IRP_ENDIF + + if (exists) then + if (mpi_master) then + call ezfio_get_mo_basis_mo_coef_aux(mo_coef_aux) + write(*,*) 'Read mo_coef_aux' + endif + IRP_IF MPI + call MPI_BCAST(mo_coef_aux, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_aux with MPI' + endif + IRP_ENDIF + else + ! Orthonormalized AO basis + do i = 1, mo_num + do j = 1, ao_num + mo_coef_aux(j,i) = ao_ortho_canonical_coef(j,i) + enddo + enddo + endif + +END_PROVIDER + diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index 0adc49fb..d1a044cb 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -408,24 +408,71 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) implicit none double precision, intent(in) :: r1(3), r2(3) double precision, intent(out) :: mu_val, mu_der(3) + double precision :: r(3) + double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) + double precision :: dm_tot, tmp1, tmp2, tmp3 if(j1b_type .eq. 200) then - double precision :: r(3), dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) - double precision :: dm_tot, tmp + ! + ! r = 0.5 (r1 + r2) + ! + ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho) + ! + ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx + ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) + ! + + PROVIDE mu_r_ct mu_erf - PROVIDE mu_r_ct r(1) = 0.5d0 * (r1(1) + r2(1)) r(2) = 0.5d0 * (r1(2) + r2(2)) r(3) = 0.5d0 * (r1(3) + r2(3)) call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) - dm_tot = dm_a(1) + dm_b(1) - mu_val = mu_r_ct * dsqrt(dm_tot) - tmp = 0.25d0 * mu_r_ct / dm_tot - mu_der(1) = tmp * (grad_dm_a(1,1) + grad_dm_b(1,1)) - mu_der(2) = tmp * (grad_dm_a(2,1) + grad_dm_b(2,1)) - mu_der(3) = tmp * (grad_dm_a(3,1) + grad_dm_b(3,1)) + + dm_tot = dm_a(1) + dm_b(1) + tmp1 = dsqrt(dm_tot) + tmp2 = mu_erf * dexp(-dm_tot) + + mu_val = mu_r_ct * tmp1 + tmp2 + + mu_der = 0.d0 + if(dm_tot .lt. 1d-7) return + + tmp3 = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2 + mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) + mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) + mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) + + elseif(j1b_type .eq. 201) then + + ! + ! r = 0.5 (r1 + r2) + ! + ! mu[rho(r)] = alpha rho + mu0 exp(-rho) + ! + ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx + ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) + ! + + PROVIDE mu_r_ct mu_erf + + r(1) = 0.5d0 * (r1(1) + r2(1)) + r(2) = 0.5d0 * (r1(2) + r2(2)) + r(3) = 0.5d0 * (r1(3) + r2(3)) + + call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) + + dm_tot = dm_a(1) + dm_b(1) + tmp2 = mu_erf * dexp(-dm_tot) + + mu_val = mu_r_ct * dm_tot + tmp2 + + tmp3 = 0.5d0 * (mu_r_ct - tmp2) + mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) + mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) + mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) else diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f index e5f123a7..980d12de 100644 --- a/src/tc_bi_ortho/print_tc_energy.irp.f +++ b/src/tc_bi_ortho/print_tc_energy.irp.f @@ -9,6 +9,10 @@ program print_tc_energy my_n_pt_a_grid = 50 read_wf = .True. touch read_wf + + PROVIDE j1b_type + print*, 'j1b_type = ', j1b_type + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call write_tc_energy end diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 2e14488e..fad6f1c2 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -176,7 +176,7 @@ default: 1.e-7 type: logical doc: If |true|, the integrals of the three-body jastrow are computed with cycles interface: ezfio,provider,ocaml -default: True +default: Flase [thresh_biorthog_diag] type: Threshold From 402a6e89885d4501d3101a65d0afc99590e7ce47 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sun, 7 May 2023 12:44:59 +0200 Subject: [PATCH 25/28] added jast_type 0 --- src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 68 ++++++++++++++++++++ src/ao_tc_eff_map/compute_ints_eff_pot.irp.f | 14 ++-- src/bi_ort_ints/one_e_bi_ort.irp.f | 26 ++++---- src/non_h_ints_mu/jast_deriv.irp.f | 49 ++++++++++++-- src/non_h_ints_mu/tc_integ.irp.f | 51 ++++++++++++++- src/non_h_ints_mu/total_tc_int.irp.f | 16 +++++ src/tc_keywords/EZFIO.cfg | 2 +- 7 files changed, 196 insertions(+), 30 deletions(-) diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 8196614f..722ff2ff 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -1,4 +1,72 @@ + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2 + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: tmp + double precision :: wall0, wall1 + + double precision, external :: overlap_gauss_r12_ao + + print*, ' providing int2_grad1u2_grad2u2 ...' + call wall_time(wall0) + + provide mu_erf final_grid_points j1b_pen + + int2_grad1u2_grad2u2 = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) + !$OMP DO + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = i, ao_num + + tmp = 0.d0 + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = coef_gauss_1_erf_x_2(i_fit) + + tmp += -0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j) + enddo + + int2_grad1u2_grad2u2(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + int2_grad1u2_grad2u2(j,i,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0 + +END_PROVIDER + ! --- BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)] diff --git a/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f b/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f index 7a567979..963a49a6 100644 --- a/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f +++ b/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f @@ -53,13 +53,13 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va integral_erf = ao_two_e_integral_erf(i, k, j, l) integral = integral_erf + integral_pot - if( j1b_type .eq. 1 ) then - !print *, ' j1b type 1 is added' - integral = integral + j1b_gauss_2e_j1(i, k, j, l) - elseif( j1b_type .eq. 2 ) then - !print *, ' j1b type 2 is added' - integral = integral + j1b_gauss_2e_j2(i, k, j, l) - endif + !if( j1b_type .eq. 1 ) then + ! !print *, ' j1b type 1 is added' + ! integral = integral + j1b_gauss_2e_j1(i, k, j, l) + !elseif( j1b_type .eq. 2 ) then + ! !print *, ' j1b type 2 is added' + ! integral = integral + j1b_gauss_2e_j2(i, k, j, l) + !endif if(abs(integral) < thr) then cycle diff --git a/src/bi_ort_ints/one_e_bi_ort.irp.f b/src/bi_ort_ints/one_e_bi_ort.irp.f index b0455570..5f2795f1 100644 --- a/src/bi_ort_ints/one_e_bi_ort.irp.f +++ b/src/bi_ort_ints/one_e_bi_ort.irp.f @@ -8,22 +8,22 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)] ao_one_e_integrals_tc_tot = ao_one_e_integrals - provide j1b_type + !provide j1b_type - if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then - - print *, ' do things properly !' - stop + !if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then + ! + ! print *, ' do things properly !' + ! stop - !do i = 1, ao_num - ! do j = 1, ao_num - ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & - ! + j1b_gauss_hermII (j,i) & - ! + j1b_gauss_nonherm(j,i) ) - ! enddo - !enddo + ! !do i = 1, ao_num + ! ! do j = 1, ao_num + ! ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & + ! ! + j1b_gauss_hermII (j,i) & + ! ! + j1b_gauss_nonherm(j,i) ) + ! ! enddo + ! !enddo - endif + !endif END_PROVIDER diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index d1a044cb..4cfd13d2 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -19,8 +19,12 @@ END_DOC implicit none - integer :: ipoint, jpoint - double precision :: r1(3), r2(3) + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: v1b_r1, v1b_r2, u2b_r12 + double precision :: grad1_v1b(3), grad1_u2b(3) + double precision :: dx, dy, dz + double precision, external :: j12_mu, j1b_nucl PROVIDE j1b_type PROVIDE final_grid_points_extra @@ -28,12 +32,43 @@ grad1_u12_num = 0.d0 grad1_u12_squared_num = 0.d0 - if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then + if(j1b_type .eq. 100) then - double precision :: v1b_r1, v1b_r2, u2b_r12 - double precision :: grad1_v1b(3), grad1_u2b(3) - double precision :: dx, dy, dz - double precision, external :: j12_mu, j1b_nucl + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) & + !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, & + !$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + call grad1_j12_mu(r1, r2, grad1_u2b) + + dx = grad1_u2b(1) + dy = grad1_u2b(2) + dz = grad1_u2b(3) + + grad1_u12_num(jpoint,ipoint,1) = dx + grad1_u12_num(jpoint,ipoint,2) = dy + grad1_u12_num(jpoint,ipoint,3) = dz + + grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then !$OMP PARALLEL & !$OMP DEFAULT (NONE) & diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index 784947b4..1a86a2e7 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -41,7 +41,34 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f else - if(j1b_type .eq. 3) then + if(j1b_type .eq. 0) then + + PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu + + int2_grad1_u12_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1) & + !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points & + !$OMP , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u12_ao) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) + int2_grad1_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) + int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) + int2_grad1_u12_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif(j1b_type .eq. 3) then PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b @@ -172,7 +199,27 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p PROVIDE j1b_type - if(j1b_type .eq. 3) then + if(j1b_type .eq. 0) then + + PROVIDE int2_grad1u2_grad2u2 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, int2_grad1u2_grad2u2) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif(j1b_type .eq. 3) then PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index a60c99da..450bbef0 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -11,6 +11,13 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, call wall_time(wall0) if(test_cycle_tc) then + + PROVIDE j1b_type + if(j1b_type .ne. 3) then + print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type + stop + endif + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num @@ -20,7 +27,9 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, enddo enddo enddo + else + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num @@ -30,6 +39,7 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, enddo enddo enddo + endif call wall_time(wall1) @@ -50,6 +60,12 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao if(test_cycle_tc) then + PROVIDE j1b_type + if(j1b_type .ne. 3) then + print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type + stop + endif + ao_tc_int_chemist = ao_tc_int_chemist_test else diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index fad6f1c2..41c07d0b 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -176,7 +176,7 @@ default: 1.e-7 type: logical doc: If |true|, the integrals of the three-body jastrow are computed with cycles interface: ezfio,provider,ocaml -default: Flase +default: False [thresh_biorthog_diag] type: Threshold From f985af03953b2d25e6706cac2412eb38cf544075 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sun, 7 May 2023 15:07:51 +0200 Subject: [PATCH 26/28] jast 4 added --- src/ao_many_one_e_ints/NEED | 1 + src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 62 ++-- .../grad_lapl_jmu_modif.irp.f | 62 ++-- src/ao_many_one_e_ints/listj1b.irp.f | 296 +++++++++++++----- src/non_h_ints_mu/grad_squared.irp.f | 2 +- src/non_h_ints_mu/j12_nucl_utils.irp.f | 173 +++++++--- src/non_h_ints_mu/tc_integ.irp.f | 4 +- 7 files changed, 404 insertions(+), 196 deletions(-) diff --git a/src/ao_many_one_e_ints/NEED b/src/ao_many_one_e_ints/NEED index 0d08442c..c57219cd 100644 --- a/src/ao_many_one_e_ints/NEED +++ b/src/ao_many_one_e_ints/NEED @@ -3,3 +3,4 @@ ao_two_e_ints becke_numerical_grid mo_one_e_ints dft_utils_in_r +tc_keywords diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 722ff2ff..fb4d71f3 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -25,11 +25,11 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin int2_grad1u2_grad2u2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & - !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -51,8 +51,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -94,15 +94,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n int2_grad1u2_grad2u2_j1b2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & - !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & - !$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & + !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & + !$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -121,7 +121,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += -0.25d0 * coef_fit * int_fit -! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle +! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle ! --- @@ -146,8 +146,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -188,15 +188,15 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final int2_u2_j1b2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & - !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & - !$OMP List_all_comb_b3_cent, int2_u2_j1b2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & + !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & + !$OMP List_all_comb_b3_cent, int2_u2_j1b2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -215,7 +215,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += coef_fit * int_fit -! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle +! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle ! --- @@ -240,8 +240,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index fc30cd83..25bb2f8b 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -24,12 +24,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po v_ij_erf_rk_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & - !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & - !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & + !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -51,7 +51,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) -! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle +! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle tmp += coef * (int_mu - int_coulomb) @@ -77,8 +77,8 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -112,13 +112,13 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ x_v_ij_erf_rk_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & - !$OMP tmp_x, tmp_y, tmp_z) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& - !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & - !$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & + !$OMP tmp_x, tmp_y, tmp_z) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & + !$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) -! if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle +! if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle tmp_x += coef * (ints(1) - ints_coulomb(1)) tmp_y += coef * (ints(2) - ints_coulomb(2)) @@ -175,8 +175,8 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -220,15 +220,15 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ v_ij_u_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & - !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & - !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & + !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ B_center(3) = List_all_comb_b2_cent(3,1) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) -! if(dabs(int_fit*coef) .lt. 1d-12) cycle +! if(dabs(int_fit*coef) .lt. 1d-12) cycle tmp += coef * coef_fit * int_fit @@ -280,8 +280,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index 4698cb27..c06d64bb 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -1,17 +1,34 @@ ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b2_size] +BEGIN_PROVIDER [integer, List_all_comb_b2_size] implicit none - List_all_comb_b2_size = 2**nucl_num + PROVIDE j1b_type + + if(j1b_type .eq. 3) then + + List_all_comb_b2_size = 2**nucl_num + + elseif(j1b_type .eq. 4) then + + List_all_comb_b2_size = nucl_num + 1 + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif + + print *, ' nb of linear terms in the envelope is ', List_all_comb_b2_size END_PROVIDER ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] +BEGIN_PROVIDER [integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] implicit none integer :: i, j @@ -50,57 +67,79 @@ END_PROVIDER List_all_comb_b2_expo = 0.d0 List_all_comb_b2_cent = 0.d0 - do i = 1, List_all_comb_b2_size + if(j1b_type .eq. 3) then - tmp_cent_x = 0.d0 - tmp_cent_y = 0.d0 - tmp_cent_z = 0.d0 - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - List_all_comb_b2_expo(i) += tmp_alphaj - tmp_cent_x += tmp_alphaj * nucl_coord(j,1) - tmp_cent_y += tmp_alphaj * nucl_coord(j,2) - tmp_cent_z += tmp_alphaj * nucl_coord(j,3) - enddo + do i = 1, List_all_comb_b2_size - if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle - - List_all_comb_b2_cent(1,i) = tmp_cent_x / List_all_comb_b2_expo(i) - List_all_comb_b2_cent(2,i) = tmp_cent_y / List_all_comb_b2_expo(i) - List_all_comb_b2_cent(3,i) = tmp_cent_z / List_all_comb_b2_expo(i) - enddo - - ! --- - - do i = 1, List_all_comb_b2_size - - do j = 2, nucl_num, 1 - tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - do k = 1, j-1, 1 - tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k) - - List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & - + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & - + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + tmp_cent_x = 0.d0 + tmp_cent_y = 0.d0 + tmp_cent_z = 0.d0 + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + List_all_comb_b2_expo(i) += tmp_alphaj + tmp_cent_x += tmp_alphaj * nucl_coord(j,1) + tmp_cent_y += tmp_alphaj * nucl_coord(j,2) + tmp_cent_z += tmp_alphaj * nucl_coord(j,3) enddo + + if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle + + List_all_comb_b2_cent(1,i) = tmp_cent_x / List_all_comb_b2_expo(i) + List_all_comb_b2_cent(2,i) = tmp_cent_y / List_all_comb_b2_expo(i) + List_all_comb_b2_cent(3,i) = tmp_cent_z / List_all_comb_b2_expo(i) enddo - if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle + ! --- - List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i) - enddo + do i = 1, List_all_comb_b2_size - ! --- + do j = 2, nucl_num, 1 + tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + do k = 1, j-1, 1 + tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k) - do i = 1, List_all_comb_b2_size + List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & + + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & + + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + enddo + enddo - phase = 0 - do j = 1, nucl_num - phase += List_all_comb_b2(j,i) + if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle + + List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i) enddo - List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) - enddo + ! --- + + do i = 1, List_all_comb_b2_size + + phase = 0 + do j = 1, nucl_num + phase += List_all_comb_b2(j,i) + enddo + + List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) + enddo + + elseif(j1b_type .eq. 4) then + + List_all_comb_b2_coef( 1) = 1.d0 + List_all_comb_b2_expo( 1) = 0.d0 + List_all_comb_b2_cent(1:3,1) = 0.d0 + do i = 1, nucl_num + List_all_comb_b2_coef( i+1) = -1.d0 + List_all_comb_b2_expo( i+1) = j1b_pen( i) + List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1) + List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2) + List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3) + enddo + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif !print *, ' coeff, expo & cent of list b2' !do i = 1, List_all_comb_b2_size @@ -115,14 +154,31 @@ END_PROVIDER BEGIN_PROVIDER [ integer, List_all_comb_b3_size] implicit none + double precision :: tmp - List_all_comb_b3_size = 3**nucl_num + if(j1b_type .eq. 3) then + + List_all_comb_b3_size = 3**nucl_num + + elseif(j1b_type .eq. 4) then + + tmp = 0.5d0 * dble(nucl_num) * (dble(nucl_num) + 3.d0) + List_all_comb_b3_size = int(tmp) + 1 + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif + + print *, ' nb of linear terms in the square of the envelope is ', List_all_comb_b3_size END_PROVIDER ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] +BEGIN_PROVIDER [integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] implicit none integer :: i, j, ii, jj @@ -162,7 +218,11 @@ END_PROVIDER implicit none integer :: i, j, k, phase + integer :: ii double precision :: tmp_alphaj, tmp_alphak, facto + double precision :: tmp1, tmp2, tmp3, tmp4 + double precision :: xi, yi, zi, xj, yj, zj + double precision :: dx, dy, dz, r2 provide j1b_pen @@ -170,60 +230,126 @@ END_PROVIDER List_all_comb_b3_expo = 0.d0 List_all_comb_b3_cent = 0.d0 - do i = 1, List_all_comb_b3_size + if(j1b_type .eq. 3) then - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - List_all_comb_b3_expo(i) += tmp_alphaj - List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) - List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) - List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3) + do i = 1, List_all_comb_b3_size + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) + List_all_comb_b3_expo(i) += tmp_alphaj + List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) + List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) + List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3) + + enddo + + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) + + List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i) + List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i) + List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i) enddo - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) + ! --- - List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i) - List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i) - List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i) - enddo + do i = 1, List_all_comb_b3_size - ! --- + do j = 2, nucl_num, 1 + tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) + do k = 1, j-1, 1 + tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k) - do i = 1, List_all_comb_b3_size + List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & + + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & + + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + enddo + enddo - do j = 2, nucl_num, 1 - tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - do k = 1, j-1, 1 - tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k) + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & - + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & - + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_b3_size + + facto = 1.d0 + phase = 0 + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b3(j,i)) + + facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj)) + phase += List_all_comb_b3(j,i) + enddo + + List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) + enddo + + elseif(j1b_type .eq. 4) then + + ii = 1 + List_all_comb_b3_coef( ii) = 1.d0 + List_all_comb_b3_expo( ii) = 0.d0 + List_all_comb_b3_cent(1:3,ii) = 0.d0 + + do i = 1, nucl_num + ii = ii + 1 + List_all_comb_b3_coef( ii) = -2.d0 + List_all_comb_b3_expo( ii) = j1b_pen( i) + List_all_comb_b3_cent(1,ii) = nucl_coord(i,1) + List_all_comb_b3_cent(2,ii) = nucl_coord(i,2) + List_all_comb_b3_cent(3,ii) = nucl_coord(i,3) + enddo + + do i = 1, nucl_num + ii = ii + 1 + List_all_comb_b3_coef( ii) = 1.d0 + List_all_comb_b3_expo( ii) = 2.d0 * j1b_pen(i) + List_all_comb_b3_cent(1,ii) = nucl_coord(i,1) + List_all_comb_b3_cent(2,ii) = nucl_coord(i,2) + List_all_comb_b3_cent(3,ii) = nucl_coord(i,3) + enddo + + do i = 1, nucl_num-1 + + tmp1 = j1b_pen(i) + + xi = nucl_coord(i,1) + yi = nucl_coord(i,2) + zi = nucl_coord(i,3) + + do j = i+1, nucl_num + + tmp2 = j1b_pen(j) + tmp3 = tmp1 + tmp2 + tmp4 = 1.d0 / tmp3 + + xj = nucl_coord(j,1) + yj = nucl_coord(j,2) + zj = nucl_coord(j,3) + + dx = xi - xj + dy = yi - yj + dz = zi - zj + r2 = dx*dx + dy*dy + dz*dz + + ii = ii + 1 + List_all_comb_b3_coef( ii) = dexp(-tmp1*tmp2*tmp4*r2) + List_all_comb_b3_expo( ii) = tmp3 + List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj) + List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj) + List_all_comb_b3_cent(3,ii) = tmp4 * (tmp1 * zi + tmp2 * zj) enddo enddo - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + else - List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i) - enddo + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop - ! --- - - do i = 1, List_all_comb_b3_size - - facto = 1.d0 - phase = 0 - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b3(j,i)) - - facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj)) - phase += List_all_comb_b3(j,i) - enddo - - List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) - enddo + endif !print *, ' coeff, expo & cent of list b3' !do i = 1, List_all_comb_b3_size diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index ea4bd36c..c93f62de 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -267,7 +267,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid) ] +BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid)] implicit none integer :: ipoint, i, j diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index c9f62b18..63a9db4b 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -8,79 +8,160 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)] double precision :: x, y, z, dx, dy, dz double precision :: a, d, e, fact_r - do ipoint = 1, n_points_final_grid + if(j1b_type .eq. 3) then - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)] - fact_r = 1.d0 - do j = 1, nucl_num - a = j1b_pen(j) - dx = x - nucl_coord(j,1) - dy = y - nucl_coord(j,2) - dz = z - nucl_coord(j,3) - d = dx*dx + dy*dy + dz*dz - e = 1.d0 - dexp(-a*d) + do ipoint = 1, n_points_final_grid - fact_r = fact_r * e + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 1.d0 + do j = 1, nucl_num + a = j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + d = dx*dx + dy*dy + dz*dz + e = 1.d0 - dexp(-a*d) + + fact_r = fact_r * e + enddo + + v_1b(ipoint) = fact_r enddo - v_1b(ipoint) = fact_r - enddo + elseif(j1b_type .eq. 4) then + + ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2) + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 1.d0 + do j = 1, nucl_num + a = j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + d = dx*dx + dy*dy + dz*dz + + fact_r = fact_r - dexp(-a*d) + enddo + + v_1b(ipoint) = fact_r + enddo + + else + + print*, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)] +BEGIN_PROVIDER [double precision, v_1b_grad, (3, n_points_final_grid)] implicit none integer :: ipoint, i, j, phase - double precision :: x, y, z, dx, dy, dz + double precision :: x, y, z, dx, dy, dz, r2 double precision :: a, d, e double precision :: fact_x, fact_y, fact_z double precision :: ax_der, ay_der, az_der, a_expo - do ipoint = 1, n_points_final_grid + PROVIDE j1b_type - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + if(j1b_type .eq. 3) then - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, List_all_comb_b2_size + ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)] - phase = 0 - a_expo = 0.d0 - ax_der = 0.d0 - ay_der = 0.d0 - az_der = 0.d0 + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, List_all_comb_b2_size + + phase = 0 + a_expo = 0.d0 + ax_der = 0.d0 + ay_der = 0.d0 + az_der = 0.d0 + do j = 1, nucl_num + a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + + phase += List_all_comb_b2(j,i) + a_expo += a * (dx*dx + dy*dy + dz*dz) + ax_der += a * dx + ay_der += a * dy + az_der += a * dz + enddo + e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) + + fact_x += e * ax_der + fact_y += e * ay_der + fact_z += e * az_der + enddo + + v_1b_grad(1,ipoint) = fact_x + v_1b_grad(2,ipoint) = fact_y + v_1b_grad(3,ipoint) = fact_z + enddo + + elseif(j1b_type .eq. 4) then + + ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2) + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + ax_der = 0.d0 + ay_der = 0.d0 + az_der = 0.d0 do j = 1, nucl_num - a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + dx = x - nucl_coord(j,1) dy = y - nucl_coord(j,2) dz = z - nucl_coord(j,3) - - phase += List_all_comb_b2(j,i) - a_expo += a * (dx*dx + dy*dy + dz*dz) - ax_der += a * dx - ay_der += a * dy - az_der += a * dz - enddo - e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) + r2 = dx*dx + dy*dy + dz*dz - fact_x += e * ax_der - fact_y += e * ay_der - fact_z += e * az_der + a = j1b_pen(j) + e = a * dexp(-a * r2) + + ax_der += e * dx + ay_der += e * dy + az_der += e * dz + enddo + + v_1b_grad(1,ipoint) = 2.d0 * ax_der + v_1b_grad(2,ipoint) = 2.d0 * ay_der + v_1b_grad(3,ipoint) = 2.d0 * az_der enddo - v_1b_grad(1,ipoint) = fact_x - v_1b_grad(2,ipoint) = fact_y - v_1b_grad(3,ipoint) = fact_z - enddo + else + + print*, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif END_PROVIDER diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index 1a86a2e7..17c24d4b 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -68,7 +68,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f !$OMP END DO !$OMP END PARALLEL - elseif(j1b_type .eq. 3) then + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b @@ -219,7 +219,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p !$OMP END DO !$OMP END PARALLEL - elseif(j1b_type .eq. 3) then + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 From 6d4975041201aab26f3dd6a4f019b6eeaa189822 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sun, 7 May 2023 18:52:03 +0200 Subject: [PATCH 27/28] rotate lr angles separately --- src/tc_scf/minimize_tc_angles.irp.f | 4 ++++ src/tc_scf/tc_scf.irp.f | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/tc_scf/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f index 1363e62b..5d3ff7f0 100644 --- a/src/tc_scf/minimize_tc_angles.irp.f +++ b/src/tc_scf/minimize_tc_angles.irp.f @@ -8,6 +8,10 @@ program print_angles my_n_pt_a_grid = 50 touch my_n_pt_r_grid my_n_pt_a_grid ! call sort_by_tc_fock + + ! TODO + ! check if rotations of orbitals affect the TC energy + ! and refuse the step call minimize_tc_orb_angles end diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 88ddd26c..2485ee8b 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -53,7 +53,7 @@ program tc_scf stop endif - call minimize_tc_orb_angles() + !call minimize_tc_orb_angles() endif From b9732d78dea26e5abf850370f275ab68d8ff481e Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Mon, 8 May 2023 23:31:20 +0200 Subject: [PATCH 28/28] IPP astice: OK --- src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 2 +- src/non_h_ints_mu/debug_fit.irp.f | 156 ++++++- src/non_h_ints_mu/grad_squared.irp.f | 65 ++- src/non_h_ints_mu/j12_nucl_utils.irp.f | 54 ++- src/non_h_ints_mu/jast_deriv.irp.f | 190 +++++++- src/non_h_ints_mu/tc_integ.irp.f | 9 +- src/non_h_ints_mu/test_non_h_ints.irp.f | 458 ++++++++++++++++++- src/tc_scf/tc_scf_energy.irp.f | 2 +- 8 files changed, 894 insertions(+), 42 deletions(-) diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index fb4d71f3..7c68de75 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -164,7 +164,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f index 5995bffa..547fd198 100644 --- a/src/non_h_ints_mu/debug_fit.irp.f +++ b/src/non_h_ints_mu/debug_fit.irp.f @@ -18,13 +18,13 @@ program debug_fit PROVIDE mu_erf j1b_pen !call test_j1b_nucl() - call test_grad_j1b_nucl() + !call test_grad_j1b_nucl() !call test_lapl_j1b_nucl() !call test_list_b2() - !call test_list_b3() + call test_list_b3() - call test_fit_u() + !call test_fit_u() !call test_fit_u2() !call test_fit_ugradu() @@ -236,16 +236,25 @@ subroutine test_list_b3() integer :: ipoint double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz double precision :: r(3) - double precision, external :: j1b_nucl + double precision :: grad_num(3), eps_der, eps_lap, tmp_der, tmp_lap, i0, ip, im + double precision, external :: j1b_nucl_square print*, ' test_list_b3 ...' + eps_ij = 1d-7 + + eps_der = 1d-5 + tmp_der = 0.5d0 / eps_der + + eps_lap = 1d-4 + tmp_lap = 1.d0 / (eps_lap*eps_lap) + + ! --- + PROVIDE v_1b_list_b3 - eps_ij = 1d-7 acc_tot = 0.d0 normalz = 0.d0 - do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -253,8 +262,7 @@ subroutine test_list_b3() r(3) = final_grid_points(3,ipoint) i_exc = v_1b_list_b3(ipoint) - i_tmp = j1b_nucl(r) - i_num = i_tmp * i_tmp + i_num = j1b_nucl_square(r) acc_ij = dabs(i_exc - i_num) if(acc_ij .gt. eps_ij) then print *, ' problem in list_b3 on', ipoint @@ -267,8 +275,136 @@ subroutine test_list_b3() normalz += dabs(i_num) enddo - print*, ' acc_tot = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot on val = ', acc_tot + print*, ' normalz on val = ', normalz + + ! --- + + PROVIDE v_1b_square_grad + + acc_tot = 0.d0 + normalz = 0.d0 + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b_square_grad(ipoint,1) + r(1) = r(1) + eps_der + ip = j1b_nucl_square(r) + r(1) = r(1) - 2.d0 * eps_der + im = j1b_nucl_square(r) + r(1) = r(1) + eps_der + i_num = tmp_der * (ip - im) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad_x list_b3 on', ipoint + print *, ' r = ', r + print *, ' r2 = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3) + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + i_exc = v_1b_square_grad(ipoint,2) + r(2) = r(2) + eps_der + ip = j1b_nucl_square(r) + r(2) = r(2) - 2.d0 * eps_der + im = j1b_nucl_square(r) + r(2) = r(2) + eps_der + i_num = tmp_der * (ip - im) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad_y list_b3 on', ipoint + print *, ' r = ', r + print *, ' r2 = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3) + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + i_exc = v_1b_square_grad(ipoint,3) + r(3) = r(3) + eps_der + ip = j1b_nucl_square(r) + r(3) = r(3) - 2.d0 * eps_der + im = j1b_nucl_square(r) + r(3) = r(3) + eps_der + i_num = tmp_der * (ip - im) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad_z list_b3 on', ipoint + print *, ' r = ', r + print *, ' r2 = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3) + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot on grad = ', acc_tot + print*, ' normalz on grad = ', normalz + + ! --- + + PROVIDE v_1b_square_lapl + + acc_tot = 0.d0 + normalz = 0.d0 + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + i0 = j1b_nucl_square(r) + + i_exc = v_1b_square_lapl(ipoint) + + r(1) = r(1) + eps_lap + ip = j1b_nucl_square(r) + r(1) = r(1) - 2.d0 * eps_lap + im = j1b_nucl_square(r) + r(1) = r(1) + eps_lap + i_num = tmp_lap * (ip - 2.d0 * i0 + im) + + r(2) = r(2) + eps_lap + ip = j1b_nucl_square(r) + r(2) = r(2) - 2.d0 * eps_lap + im = j1b_nucl_square(r) + r(2) = r(2) + eps_lap + i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im) + + r(3) = r(3) + eps_lap + ip = j1b_nucl_square(r) + r(3) = r(3) - 2.d0 * eps_lap + im = j1b_nucl_square(r) + r(3) = r(3) + eps_lap + i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im) + + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in lapl list_b3 on', ipoint + print *, ' r = ', r + print *, ' r2 = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3) + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot on lapl = ', acc_tot + print*, ' normalz on lapl = ', normalz + + ! --- return end subroutine test_list_b3 diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index c93f62de..b4ddb606 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi ! ! if J(r1,r2) = u12 x v1 x v2 ! - ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2 v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2 v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ] + ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2 v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)) + u12^2 v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ] ! = -0.25 x v1^2 \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2 ! + -0.50 x (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2) u12^2 v2^2 ! + -1.00 x v1 (grad_1 v1) \int r2 \phi_i(2) \phi_j(2) (grad_1 u12) v2^2 @@ -358,7 +358,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao implicit none integer :: ipoint, i, j, k, l - double precision :: weight1, ao_ik_r, ao_i_r + double precision :: weight1, ao_k_r, ao_i_r + double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq double precision :: time0, time1 double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:) @@ -373,16 +374,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao else + ! --- + PROVIDE int2_grad1_u12_square_ao allocate(b_mat(n_points_final_grid,ao_num,ao_num)) b_mat = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, k, ipoint) & - !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) - !$OMP DO SCHEDULE (static) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) do i = 1, ao_num do k = 1, ao_num do ipoint = 1, n_points_final_grid @@ -390,13 +393,57 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL tc_grad_square_ao = 0.d0 call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & , 0.d0, tc_grad_square_ao, ao_num*ao_num) + + ! --- + + if((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then + + ! an additional term is added here directly instead of + ! being added in int2_grad1_u12_square_ao for performance + ! note that the factor + + PROVIDE int2_u2_j1b2 + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, & + !$OMP v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + + weight1 = 0.25d0 * final_weight_at_r_vector(ipoint) + + ao_i_r = aos_in_r_array_transp(ipoint,i) + ao_k_r = aos_in_r_array_transp(ipoint,k) + + b_mat(ipoint,k,i) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) ) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + , 1.d0, tc_grad_square_ao, ao_num*ao_num) + endif + + ! --- + deallocate(b_mat) call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num) diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index 63a9db4b..079cb388 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -59,7 +59,7 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)] else - print*, 'j1b_type = ', j1b_pen, 'is not implemented' + print*, 'j1b_type = ', j1b_pen, 'is not implemented for v_1b' stop endif @@ -172,7 +172,7 @@ BEGIN_PROVIDER [ double precision, v_1b_lapl, (n_points_final_grid)] implicit none integer :: ipoint, i, j, phase double precision :: x, y, z, dx, dy, dz - double precision :: a, d, e, b + double precision :: a, e, b double precision :: fact_r double precision :: ax_der, ay_der, az_der, a_expo @@ -283,6 +283,56 @@ BEGIN_PROVIDER [ double precision, v_1b_list_b3, (n_points_final_grid)] END_PROVIDER +! --- + + BEGIN_PROVIDER [double precision, v_1b_square_grad, (n_points_final_grid,3)] +&BEGIN_PROVIDER [double precision, v_1b_square_lapl, (n_points_final_grid) ] + + implicit none + integer :: ipoint, i + double precision :: x, y, z, dx, dy, dz, r2 + double precision :: coef, expo, a_expo, tmp + double precision :: fact_x, fact_y, fact_z, fact_r + + PROVIDE List_all_comb_b3_coef List_all_comb_b3_expo List_all_comb_b3_cent + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + fact_r = 0.d0 + do i = 1, List_all_comb_b3_size + + coef = List_all_comb_b3_coef(i) + expo = List_all_comb_b3_expo(i) + + dx = x - List_all_comb_b3_cent(1,i) + dy = y - List_all_comb_b3_cent(2,i) + dz = z - List_all_comb_b3_cent(3,i) + r2 = dx * dx + dy * dy + dz * dz + + a_expo = expo * r2 + tmp = coef * expo * dexp(-a_expo) + + fact_x += tmp * dx + fact_y += tmp * dy + fact_z += tmp * dz + fact_r += tmp * (3.d0 - 2.d0 * a_expo) + enddo + + v_1b_square_grad(ipoint,1) = -2.d0 * fact_x + v_1b_square_grad(ipoint,2) = -2.d0 * fact_y + v_1b_square_grad(ipoint,3) = -2.d0 * fact_z + v_1b_square_lapl(ipoint) = -2.d0 * fact_r + enddo + +END_PROVIDER + ! --- double precision function j12_mu_r12(r12) diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index 4cfd13d2..cbd0b406 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -164,7 +164,7 @@ double precision function j12_mu(r1, r2) double precision, intent(in) :: r1(3), r2(3) double precision :: mu_tmp, r12 - if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + (r1(2) - r2(2)) * (r1(2) - r2(2)) & @@ -175,7 +175,7 @@ double precision function j12_mu(r1, r2) else - print *, ' j1b_type = ', j1b_type, 'not implemented yet' + print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu' stop endif @@ -196,7 +196,7 @@ subroutine grad1_j12_mu(r1, r2, grad) grad = 0.d0 - if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then dx = r1(1) - r2(1) dy = r1(2) - r2(2) @@ -252,7 +252,7 @@ double precision function j1b_nucl(r) integer :: i double precision :: a, d, e, x, y, z - if(j1b_type .eq. 102) then + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then j1b_nucl = 1.d0 do i = 1, nucl_num @@ -263,7 +263,7 @@ double precision function j1b_nucl(r) j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d)) enddo - elseif(j1b_type .eq. 103) then + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then j1b_nucl = 1.d0 do i = 1, nucl_num @@ -275,7 +275,7 @@ double precision function j1b_nucl(r) j1b_nucl = j1b_nucl * e enddo - elseif(j1b_type .eq. 104) then + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then j1b_nucl = 1.d0 do i = 1, nucl_num @@ -286,7 +286,7 @@ double precision function j1b_nucl(r) j1b_nucl = j1b_nucl - dexp(-a*d) enddo - elseif(j1b_type .eq. 105) then + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then j1b_nucl = 1.d0 do i = 1, nucl_num @@ -300,7 +300,7 @@ double precision function j1b_nucl(r) else - print *, ' j1b_type = ', j1b_type, 'not implemented yet' + print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl' stop endif @@ -310,6 +310,75 @@ end function j1b_nucl ! --- +double precision function j1b_nucl_square(r) + + implicit none + double precision, intent(in) :: r(3) + integer :: i + double precision :: a, d, e, x, y, z + + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl_square = j1b_nucl_square - dexp(-a*dsqrt(d)) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + e = 1.d0 - dexp(-a*d) + j1b_nucl_square = j1b_nucl_square * e + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl_square = j1b_nucl_square - dexp(-a*d) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + j1b_nucl_square = j1b_nucl_square - dexp(-a*d*d) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_square' + stop + + endif + + return +end function j1b_nucl_square + +! --- + subroutine grad1_j1b_nucl(r, grad) implicit none @@ -321,7 +390,7 @@ subroutine grad1_j1b_nucl(r, grad) double precision :: fact_x, fact_y, fact_z double precision :: ax_der, ay_der, az_der, a_expo - if(j1b_type .eq. 102) then + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then fact_x = 0.d0 fact_y = 0.d0 @@ -343,7 +412,7 @@ subroutine grad1_j1b_nucl(r, grad) grad(2) = fact_y grad(3) = fact_z - elseif(j1b_type .eq. 103) then + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then x = r(1) y = r(2) @@ -382,7 +451,7 @@ subroutine grad1_j1b_nucl(r, grad) grad(2) = fact_y grad(3) = fact_z - elseif(j1b_type .eq. 104) then + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then fact_x = 0.d0 fact_y = 0.d0 @@ -404,7 +473,7 @@ subroutine grad1_j1b_nucl(r, grad) grad(2) = 2.d0 * fact_y grad(3) = 2.d0 * fact_z - elseif(j1b_type .eq. 105) then + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then fact_x = 0.d0 fact_y = 0.d0 @@ -428,7 +497,7 @@ subroutine grad1_j1b_nucl(r, grad) else - print *, ' j1b_type = ', j1b_type, 'not implemented yet' + print *, ' j1b_type = ', j1b_type, 'not implemented for grad1_j1b_nucl' stop endif @@ -520,3 +589,98 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) end subroutine mu_r_val_and_grad ! --- + +subroutine grad1_j1b_nucl_square_num(r1, grad) + + implicit none + double precision, intent(in) :: r1(3) + double precision, intent(out) :: grad(3) + double precision :: r(3), eps, tmp_eps, vp, vm + double precision, external :: j1b_nucl_square + + eps = 1d-5 + tmp_eps = 0.5d0 / eps + + r(1:3) = r1(1:3) + + r(1) = r(1) + eps + vp = j1b_nucl_square(r) + r(1) = r(1) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(1) = r(1) + eps + grad(1) = tmp_eps * (vp - vm) + + r(2) = r(2) + eps + vp = j1b_nucl_square(r) + r(2) = r(2) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(2) = r(2) + eps + grad(2) = tmp_eps * (vp - vm) + + r(3) = r(3) + eps + vp = j1b_nucl_square(r) + r(3) = r(3) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(3) = r(3) + eps + grad(3) = tmp_eps * (vp - vm) + + return +end subroutine grad1_j1b_nucl_square_num + +! --- + +subroutine grad1_j12_mu_square_num(r1, r2, grad) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: grad(3) + double precision :: r(3) + double precision :: eps, tmp_eps, vp, vm + double precision, external :: j12_mu_square + + eps = 1d-5 + tmp_eps = 0.5d0 / eps + + r(1:3) = r1(1:3) + + r(1) = r(1) + eps + vp = j12_mu_square(r, r2) + r(1) = r(1) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(1) = r(1) + eps + grad(1) = tmp_eps * (vp - vm) + + r(2) = r(2) + eps + vp = j12_mu_square(r, r2) + r(2) = r(2) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(2) = r(2) + eps + grad(2) = tmp_eps * (vp - vm) + + r(3) = r(3) + eps + vp = j12_mu_square(r, r2) + r(3) = r(3) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(3) = r(3) + eps + grad(3) = tmp_eps * (vp - vm) + + return +end subroutine grad1_j12_mu_square_num + +! --- + +double precision function j12_mu_square(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, external :: j12_mu + + j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2) + + return +end function j12_mu_square + +! --- + diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index 17c24d4b..efa9be43 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -221,18 +221,21 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then - PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 + ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance + !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 + PROVIDE u12sq_j1bsq grad12_j12 int2_grad1_u12_square_ao = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, j, ipoint) & - !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12) !$OMP DO SCHEDULE (static) do ipoint = 1, n_points_final_grid do j = 1, ao_num do i = 1, ao_num - int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) + !int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) + int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) enddo enddo enddo diff --git a/src/non_h_ints_mu/test_non_h_ints.irp.f b/src/non_h_ints_mu/test_non_h_ints.irp.f index c535d0c5..a6e0a311 100644 --- a/src/non_h_ints_mu/test_non_h_ints.irp.f +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -1,15 +1,25 @@ program test_non_h - implicit none + + implicit none + my_grid_becke = .True. my_n_pt_r_grid = 50 my_n_pt_a_grid = 74 + !my_n_pt_r_grid = 400 + !my_n_pt_a_grid = 974 + ! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_a_grid = 26 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid -!call routine_grad_squared - call routine_fit + + !call routine_grad_squared + !call routine_fit + + call test_ipp() end +! --- + subroutine routine_lapl_grad implicit none integer :: i,j,k,l @@ -100,3 +110,445 @@ subroutine routine_fit enddo end + + +subroutine test_ipp() + + implicit none + integer :: i, j, k, l, ipoint + double precision :: accu, norm, diff, old, new, eps, int_num + double precision :: weight1, ao_i_r, ao_k_r + double precision, allocatable :: b_mat(:,:,:), I1(:,:,:,:), I2(:,:,:,:) + + eps = 1d-7 + + allocate(b_mat(n_points_final_grid,ao_num,ao_num)) + b_mat = 0.d0 + + ! --- + + ! first way + + allocate(I1(ao_num,ao_num,ao_num,ao_num)) + I1 = 0.d0 + + PROVIDE u12_grad1_u12_j1b_grad1_j1b + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , u12_grad1_u12_j1b_grad1_j1b(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + , 0.d0, I1, ao_num*ao_num) + + ! --- + + ! 2nd way + + allocate(I2(ao_num,ao_num,ao_num,ao_num)) + I2 = 0.d0 + + PROVIDE int2_u2_j1b2 + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, & + !$OMP v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + + weight1 = 0.25d0 * final_weight_at_r_vector(ipoint) + + ao_i_r = aos_in_r_array_transp(ipoint,i) + ao_k_r = aos_in_r_array_transp(ipoint,k) + + b_mat(ipoint,k,i) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) & + + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) ) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + , 0.d0, I2, ao_num*ao_num) + + ! --- + + deallocate(b_mat) + + accu = 0.d0 + norm = 0.d0 + do i = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + do j = 1, ao_num + + old = I1(j,l,k,i) + new = I2(j,l,k,i) + + !print*, l, k, j, i + !print*, old, new + + diff = new - old + if(dabs(diff) .gt. eps) then + print*, ' problem on :', j, l, k, i + print*, ' diff = ', diff + print*, ' old value = ', old + print*, ' new value = ', new + call I_grade_gradu_naive1(i, j, k, l, int_num) + print*, ' full num1 = ', int_num + call I_grade_gradu_naive2(i, j, k, l, int_num) + print*, ' full num2 = ', int_num + call I_grade_gradu_naive3(i, j, k, l, int_num) + print*, ' full num3 = ', int_num + call I_grade_gradu_naive4(i, j, k, l, int_num) + print*, ' full num4 = ', int_num + call I_grade_gradu_seminaive(i, j, k, l, int_num) + print*, ' semi num = ', int_num + endif + + accu += dabs(diff) + norm += dabs(old) + enddo + enddo + enddo + enddo + + deallocate(I1, I2) + + print*, ' accu = ', accu + print*, ' norm = ', norm + + return +end subroutine test_ipp + +! --- + +subroutine I_grade_gradu_naive1(i, j, k, l, int) + + implicit none + integer, intent(in) :: i, j, k, l + double precision, intent(out) :: int + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: weight1_x, weight1_y, weight1_z + double precision :: weight2_x, weight2_y, weight2_z + double precision :: aor_i, aor_j, aor_k, aor_l + double precision :: e1_val, e2_val, e1_der(3), u12_val, u12_der(3) + double precision, external :: j1b_nucl, j12_mu + + int = 0.d0 + + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + aor_i = aos_in_r_array_transp(ipoint,i) + aor_k = aos_in_r_array_transp(ipoint,k) + + e1_val = j1b_nucl(r1) + call grad1_j1b_nucl(r1, e1_der) + + weight1_x = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(1) + weight1_y = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(2) + weight1_z = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(3) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + aor_j = aos_in_r_array_extra_transp(jpoint,j) + aor_l = aos_in_r_array_extra_transp(jpoint,l) + + e2_val = j1b_nucl(r2) + + u12_val = j12_mu(r1, r2) + call grad1_j12_mu(r1, r2, u12_der) + + weight2_x = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(1) + weight2_y = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(2) + weight2_z = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(3) + + int = int - (weight1_x * weight2_x + weight1_y * weight2_y + weight1_z * weight2_z) + enddo + enddo + + return +end subroutine I_grade_gradu_naive1 + +! --- + +subroutine I_grade_gradu_naive2(i, j, k, l, int) + + implicit none + integer, intent(in) :: i, j, k, l + double precision, intent(out) :: int + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: weight1_x, weight1_y, weight1_z + double precision :: weight2_x, weight2_y, weight2_z + double precision :: aor_i, aor_j, aor_k, aor_l + double precision :: e1_square_der(3), e2_val, u12_square_der(3) + double precision, external :: j1b_nucl + + int = 0.d0 + + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + aor_i = aos_in_r_array_transp(ipoint,i) + aor_k = aos_in_r_array_transp(ipoint,k) + + call grad1_j1b_nucl_square_num(r1, e1_square_der) + + weight1_x = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(1) + weight1_y = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(2) + weight1_z = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(3) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + aor_j = aos_in_r_array_extra_transp(jpoint,j) + aor_l = aos_in_r_array_extra_transp(jpoint,l) + + e2_val = j1b_nucl(r2) + call grad1_j12_mu_square_num(r1, r2, u12_square_der) + + weight2_x = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(1) + weight2_y = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(2) + weight2_z = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(3) + + int = int - 0.25d0 * (weight1_x * weight2_x + weight1_y * weight2_y + weight1_z * weight2_z) + enddo + enddo + + return +end subroutine I_grade_gradu_naive2 + +! --- + +subroutine I_grade_gradu_naive3(i, j, k, l, int) + + implicit none + integer, intent(in) :: i, j, k, l + double precision, intent(out) :: int + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: weight1, weight2 + double precision :: aor_j, aor_l + double precision :: grad(3), e2_val, u12_val + double precision, external :: j1b_nucl, j12_mu + + int = 0.d0 + + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + call grad1_aos_ik_grad1_esquare(i, k, r1, grad) + + weight1 = final_weight_at_r_vector(ipoint) * (grad(1) + grad(2) + grad(3)) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + aor_j = aos_in_r_array_extra_transp(jpoint,j) + aor_l = aos_in_r_array_extra_transp(jpoint,l) + + e2_val = j1b_nucl(r2) + u12_val = j12_mu(r1, r2) + + weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint) + + int = int + 0.25d0 * weight1 * weight2 + enddo + enddo + + return +end subroutine I_grade_gradu_naive3 + +! --- + +subroutine I_grade_gradu_naive4(i, j, k, l, int) + + implicit none + integer, intent(in) :: i, j, k, l + double precision, intent(out) :: int + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: weight1, weight2 + double precision :: aor_j, aor_l, aor_k, aor_i + double precision :: grad(3), e2_val, u12_val + double precision, external :: j1b_nucl, j12_mu + + int = 0.d0 + + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + aor_i = aos_in_r_array_transp(ipoint,i) + aor_k = aos_in_r_array_transp(ipoint,k) + + weight1 = final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) ) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + aor_j = aos_in_r_array_extra_transp(jpoint,j) + aor_l = aos_in_r_array_extra_transp(jpoint,l) + + e2_val = j1b_nucl(r2) + u12_val = j12_mu(r1, r2) + + weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint) + + int = int + 0.25d0 * weight1 * weight2 + enddo + enddo + + return +end subroutine I_grade_gradu_naive4 + +! --- + +subroutine I_grade_gradu_seminaive(i, j, k, l, int) + + implicit none + integer, intent(in) :: i, j, k, l + double precision, intent(out) :: int + integer :: ipoint + double precision :: r1(3) + double precision :: weight1 + double precision :: aor_i, aor_k + + int = 0.d0 + + do ipoint = 1, n_points_final_grid ! r1 + + aor_i = aos_in_r_array_transp(ipoint,i) + aor_k = aos_in_r_array_transp(ipoint,k) + + weight1 = 0.25d0 * final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) & + + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) ) + + int = int + weight1 * int2_u2_j1b2(j,l,ipoint) + enddo + + return +end subroutine I_grade_gradu_seminaive + +! --- + +subroutine aos_ik_grad1_esquare(i, k, r1, val) + + implicit none + integer, intent(in) :: i, k + double precision, intent(in) :: r1(3) + double precision, intent(out) :: val(3) + double precision :: tmp + double precision :: der(3), aos_array(ao_num), aos_grad_array(3,ao_num) + + call give_all_aos_and_grad_at_r(r1, aos_array, aos_grad_array) + call grad1_j1b_nucl_square_num(r1, der) + + tmp = aos_array(i) * aos_array(k) + val(1) = tmp * der(1) + val(2) = tmp * der(2) + val(3) = tmp * der(3) + + return +end subroutine phi_ik_grad1_esquare + +! --- + +subroutine grad1_aos_ik_grad1_esquare(i, k, r1, grad) + + implicit none + integer, intent(in) :: i, k + double precision, intent(in) :: r1(3) + double precision, intent(out) :: grad(3) + double precision :: r(3), eps, tmp_eps, val_p(3), val_m(3) + + eps = 1d-5 + tmp_eps = 0.5d0 / eps + + r(1:3) = r1(1:3) + + r(1) = r(1) + eps + call aos_ik_grad1_esquare(i, k, r, val_p) + r(1) = r(1) - 2.d0 * eps + call aos_ik_grad1_esquare(i, k, r, val_m) + r(1) = r(1) + eps + grad(1) = tmp_eps * (val_p(1) - val_m(1)) + + r(2) = r(2) + eps + call aos_ik_grad1_esquare(i, k, r, val_p) + r(2) = r(2) - 2.d0 * eps + call aos_ik_grad1_esquare(i, k, r, val_m) + r(2) = r(2) + eps + grad(2) = tmp_eps * (val_p(2) - val_m(2)) + + r(3) = r(3) + eps + call aos_ik_grad1_esquare(i, k, r, val_p) + r(3) = r(3) - 2.d0 * eps + call aos_ik_grad1_esquare(i, k, r, val_m) + r(3) = r(3) + eps + grad(3) = tmp_eps * (val_p(3) - val_m(3)) + + return +end subroutine grad1_aos_ik_grad1_esquare + +! --- + + + + + + + diff --git a/src/tc_scf/tc_scf_energy.irp.f b/src/tc_scf/tc_scf_energy.irp.f index 5c643f19..833b48aa 100644 --- a/src/tc_scf/tc_scf_energy.irp.f +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ double precision, TC_HF_energy] + BEGIN_PROVIDER [ double precision, TC_HF_energy ] &BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy] &BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy]