From de288449f58a54893cf1647faa8b00116303e7bc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Apr 2024 13:45:51 +0200 Subject: [PATCH 1/6] Fix dos files in qp_create --- ocaml/Atom.ml | 13 +++++++++---- ocaml/Molecule.ml | 12 +++++++++++- ocaml/Point3d.ml | 4 +++- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index d02b20d8..49e788e8 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -22,10 +22,15 @@ let of_string ~units s = } | [ name; x; y; z ] -> let e = Element.of_string name in - { element = e ; - charge = Element.to_charge e; - coord = Point3d.of_string ~units (String.concat " " [x; y; z]) - } + begin + try + { element = e ; + charge = Element.to_charge e; + coord = Point3d.of_string ~units (String.concat " " [x; y; z]) + } + with + | err -> (Printf.eprintf "name = \"%s\"\nxyz = (%s,%s,%s)\n%!" name x y z ; raise err) + end | _ -> raise (AtomError s) diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 603244c8..3771b6f9 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -142,13 +142,21 @@ let of_xyz_string result +let regexp_r = Str.regexp {| |} +let regexp_t = Str.regexp {| |} + let of_xyz_file ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) ?(units=Units.Angstrom) filename = let lines = - match Io_ext.input_lines filename with + Io_ext.input_lines filename + |> List.map (fun s -> Str.global_replace regexp_r "" s) + |> List.map (fun s -> Str.global_replace regexp_t " " s) + in + let lines = + match lines with | natoms :: title :: rest -> let natoms = try @@ -173,6 +181,8 @@ let of_zmt_file ?(units=Units.Angstrom) filename = Io_ext.read_all filename + |> Str.global_replace regexp_r "" + |> Str.global_replace regexp_t " " |> Zmatrix.of_string |> Zmatrix.to_xyz_string |> of_xyz_string ~charge ~multiplicity ~units diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 57b02bfe..4df375e4 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -24,7 +24,9 @@ let of_string ~units s = let l = s |> String_ext.split ~on:' ' |> List.filter (fun x -> x <> "") - |> list_map float_of_string + |> list_map (fun x -> + try float_of_string x with + | Failure msg -> (Printf.eprintf "Bad string: \"%s\"\n%!" x ; failwith msg) ) |> Array.of_list in { x = l.(0) *. f ; From 1c2b737586eba60cfec15ce8c452bdff727c70b9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 May 2024 16:05:13 +0200 Subject: [PATCH 2/6] Fixed Warning with nproc --- src/utils/util.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index 97cbde67..de01656b 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -327,12 +327,12 @@ subroutine wall_time(t) end BEGIN_PROVIDER [ integer, nproc ] + use omp_lib implicit none BEGIN_DOC ! Number of current OpenMP threads END_DOC - integer, external :: omp_get_num_threads nproc = 1 !$OMP PARALLEL !$OMP MASTER From 425e7e4ee0ac740220bb921ba7a607836b1acffe Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 May 2024 16:22:01 +0200 Subject: [PATCH 3/6] Changed symetric_fock_tc into symmetric_fock_tc --- plugins/local/tc_keywords/EZFIO.cfg | 2 +- plugins/local/tc_scf/fock_hermit.irp.f | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg index bc691fc3..e0776136 100644 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ b/plugins/local/tc_keywords/EZFIO.cfg @@ -106,7 +106,7 @@ doc: If |true|, the MO basis is assumed to be bi-orthonormal interface: ezfio,provider,ocaml default: True -[symetric_fock_tc] +[symmetric_fock_tc] type: logical doc: If |true|, using F+F^t as Fock TC interface: ezfio,provider,ocaml diff --git a/plugins/local/tc_scf/fock_hermit.irp.f b/plugins/local/tc_scf/fock_hermit.irp.f index 5a51b324..3460157e 100644 --- a/plugins/local/tc_scf/fock_hermit.irp.f +++ b/plugins/local/tc_scf/fock_hermit.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, good_hermit_tc_fock_mat, (mo_num, mo_num)] BEGIN_DOC -! good_hermit_tc_fock_mat = Hermitian Upper triangular Fock matrix +! good_hermit_tc_fock_mat = Hermitian Upper triangular Fock matrix ! ! The converged eigenvectors of such matrix yield to orthonormal vectors satisfying the left Brillouin theorem END_DOC @@ -14,11 +14,11 @@ BEGIN_PROVIDER [ double precision, good_hermit_tc_fock_mat, (mo_num, mo_num)] good_hermit_tc_fock_mat = Fock_matrix_tc_mo_tot do j = 1, mo_num do i = 1, j-1 - good_hermit_tc_fock_mat(i,j) = Fock_matrix_tc_mo_tot(j,i) + good_hermit_tc_fock_mat(i,j) = Fock_matrix_tc_mo_tot(j,i) enddo enddo -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [ double precision, hermit_average_tc_fock_mat, (mo_num, mo_num)] @@ -35,7 +35,7 @@ BEGIN_PROVIDER [ double precision, hermit_average_tc_fock_mat, (mo_num, mo_num)] enddo enddo -END_PROVIDER +END_PROVIDER ! --- @@ -44,13 +44,13 @@ BEGIN_PROVIDER [ double precision, grad_hermit] BEGIN_DOC ! square of gradient of the energy END_DOC - if(symetric_fock_tc)then + if(symmetric_fock_tc)then grad_hermit = grad_hermit_average_tc_fock_mat else grad_hermit = grad_good_hermit_tc_fock_mat endif -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [ double precision, grad_good_hermit_tc_fock_mat] implicit none @@ -64,7 +64,7 @@ BEGIN_PROVIDER [ double precision, grad_good_hermit_tc_fock_mat] grad_good_hermit_tc_fock_mat += dabs(good_hermit_tc_fock_mat(i,j)) enddo enddo -END_PROVIDER +END_PROVIDER ! --- @@ -80,7 +80,7 @@ BEGIN_PROVIDER [ double precision, grad_hermit_average_tc_fock_mat] grad_hermit_average_tc_fock_mat += dabs(hermit_average_tc_fock_mat(i,j)) enddo enddo -END_PROVIDER +END_PROVIDER ! --- @@ -95,8 +95,8 @@ subroutine save_good_hermit_tc_eigvectors() sign = 1 label = "Canonical" output = .False. - - if(symetric_fock_tc)then + + if(symmetric_fock_tc)then call mo_as_eigvectors_of_mo_matrix(hermit_average_tc_fock_mat, mo_num, mo_num, label, sign, output) else call mo_as_eigvectors_of_mo_matrix(good_hermit_tc_fock_mat, mo_num, mo_num, label, sign, output) From 13785b267c36319925ffa72ebe42399fa932ffae Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 3 May 2024 11:34:30 +0200 Subject: [PATCH 4/6] fixed a bug in src/scf_utils/roothaan_hall_scf.irp.f --- .../extra_grid_vector.irp.f | 20 +++++++++---------- .../grid_becke_vector.irp.f | 20 +++++++++---------- src/scf_utils/roothaan_hall_scf.irp.f | 2 +- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/becke_numerical_grid/extra_grid_vector.irp.f b/src/becke_numerical_grid/extra_grid_vector.irp.f index 16a52dc6..44fc4435 100644 --- a/src/becke_numerical_grid/extra_grid_vector.irp.f +++ b/src/becke_numerical_grid/extra_grid_vector.irp.f @@ -71,16 +71,16 @@ END_PROVIDER index_final_points_extra(3,i_count) = j index_final_points_extra_reverse(k,i,j) = i_count - if(final_weight_at_r_vector_extra(i_count) .lt. 0.d0) then - print *, ' !!! WARNING !!!' - print *, ' negative weight !!!!' - print *, i_count, final_weight_at_r_vector_extra(i_count) - if(dabs(final_weight_at_r_vector_extra(i_count)) .lt. 1d-10) then - final_weight_at_r_vector_extra(i_count) = 0.d0 - else - stop - endif - endif +! if(final_weight_at_r_vector_extra(i_count) .lt. 0.d0) then +! print *, ' !!! WARNING !!!' +! print *, ' negative weight !!!!' +! print *, i_count, final_weight_at_r_vector_extra(i_count) +! if(dabs(final_weight_at_r_vector_extra(i_count)) .lt. 1d-10) then +! final_weight_at_r_vector_extra(i_count) = 0.d0 +! else +! stop +! endif +! endif enddo enddo enddo diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index c35918c3..7097dbb3 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -68,16 +68,16 @@ END_PROVIDER index_final_points(3,i_count) = j index_final_points_reverse(k,i,j) = i_count - if(final_weight_at_r_vector(i_count) .lt. 0.d0) then - print *, ' !!! WARNING !!!' - print *, ' negative weight !!!!' - print *, i_count, final_weight_at_r_vector(i_count) - if(dabs(final_weight_at_r_vector(i_count)) .lt. 1d-10) then - final_weight_at_r_vector(i_count) = 0.d0 - else - stop - endif - endif +! if(final_weight_at_r_vector(i_count) .lt. 0.d0) then +! print *, ' !!! WARNING !!!' +! print *, ' negative weight !!!!' +! print *, i_count, final_weight_at_r_vector(i_count) +! if(dabs(final_weight_at_r_vector(i_count)) .lt. 1d-10) then +! final_weight_at_r_vector(i_count) = 0.d0 +! else +! stop +! endif +! endif enddo enddo enddo diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 3f5c8549..e0fe5319 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -217,7 +217,7 @@ END_DOC do while (i Date: Mon, 6 May 2024 17:47:48 +0200 Subject: [PATCH 5/6] updated get_fci_tc_conv.sh --- scripts/get_fci_tc_conv.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/get_fci_tc_conv.sh b/scripts/get_fci_tc_conv.sh index 643f3ac0..f0c99baf 100755 --- a/scripts/get_fci_tc_conv.sh +++ b/scripts/get_fci_tc_conv.sh @@ -1,2 +1,2 @@ file=$1 -grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc +grep "Ndet,E,E+PT2,pt2_minus,pt2_plus,pt2_abs=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc From b14325fef482bdf6cb471b40edf8fa46f2aeac65 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 6 May 2024 18:21:58 +0200 Subject: [PATCH 6/6] Introducing qmckl --- plugins/local/non_h_ints_mu/NEED | 1 + plugins/local/non_h_ints_mu/deb_aos.irp.f | 49 ++++++++-- plugins/local/non_h_ints_mu/qmckl.irp.f | 104 ++++++++++++++++++++++ 3 files changed, 148 insertions(+), 6 deletions(-) diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED index 48c1c24b..5ca1d543 100644 --- a/plugins/local/non_h_ints_mu/NEED +++ b/plugins/local/non_h_ints_mu/NEED @@ -3,3 +3,4 @@ hamiltonian jastrow ao_tc_eff_map bi_ortho_mos +trexio diff --git a/plugins/local/non_h_ints_mu/deb_aos.irp.f b/plugins/local/non_h_ints_mu/deb_aos.irp.f index c9bc9c9a..86d011fb 100644 --- a/plugins/local/non_h_ints_mu/deb_aos.irp.f +++ b/plugins/local/non_h_ints_mu/deb_aos.irp.f @@ -34,21 +34,58 @@ subroutine print_aos() PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array - do ipoint = 1, n_points_final_grid - r(:) = final_grid_points(:,ipoint) - print*, r - enddo +! do ipoint = 1, n_points_final_grid +! r(:) = final_grid_points(:,ipoint) +! print*, r +! enddo +double precision :: accu_vgl(5) +double precision :: accu_vgl_nrm(5) do ipoint = 1, n_points_final_grid - r(:) = final_grid_points(:,ipoint) do i = 1, ao_num ao_val = aos_in_r_array (i,ipoint) ao_der(:) = aos_grad_in_r_array(i,ipoint,:) ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) - write(*, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + write(111, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap enddo enddo + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + ao_val = aos_in_r_array_qmckl (i,ipoint) + ao_der(:) = aos_grad_in_r_array_qmckl(i,ipoint,:) + ao_lap = aos_lapl_in_r_array_qmckl(i,ipoint) + write(222, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + enddo + enddo + + accu_vgl = 0.d0 + accu_vgl_nrm = 0.d0 + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + ao_val = aos_in_r_array (i,ipoint) + ao_der(:) = aos_grad_in_r_array(i,ipoint,:) + ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) + accu_vgl_nrm(1) += dabs(ao_val) + accu_vgl_nrm(2) += dabs(ao_der(1)) + accu_vgl_nrm(3) += dabs(ao_der(2)) + accu_vgl_nrm(4) += dabs(ao_der(3)) + accu_vgl_nrm(5) += dabs(ao_lap) + + ao_val -= aos_in_r_array_qmckl (i,ipoint) + ao_der(:) -= aos_grad_in_r_array_qmckl(i,ipoint,:) + ao_lap -= aos_lapl_in_r_array_qmckl(i,ipoint) + accu_vgl(1) += dabs(ao_val) + accu_vgl(2) += dabs(ao_der(1)) + accu_vgl(3) += dabs(ao_der(2)) + accu_vgl(4) += dabs(ao_der(3)) + accu_vgl(5) += dabs(ao_lap) + enddo + + enddo + accu_vgl(:) *= 1.d0 / accu_vgl_nrm(:) + print *, accu_vgl + return end diff --git a/plugins/local/non_h_ints_mu/qmckl.irp.f b/plugins/local/non_h_ints_mu/qmckl.irp.f index 1df80457..4d419e24 100644 --- a/plugins/local/non_h_ints_mu/qmckl.irp.f +++ b/plugins/local/non_h_ints_mu/qmckl.irp.f @@ -75,3 +75,107 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ] endif END_PROVIDER + + + BEGIN_PROVIDER [ double precision, aos_in_r_array_qmckl, (ao_num,n_points_final_grid)] +&BEGIN_PROVIDER [ double precision, aos_grad_in_r_array_qmckl, (ao_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [ double precision, aos_lapl_in_r_array_qmckl, (ao_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! AOS computed with qmckl + END_DOC + use qmckl + + integer*8 :: qmckl_ctx + integer(qmckl_exit_code) :: rc + + qmckl_ctx = qmckl_context_create() + + rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename))) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in read_trexio' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in set_electron_point' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + double precision, allocatable :: vgl(:,:,:) + allocate( vgl(ao_num,5,n_points_final_grid)) + rc = qmckl_get_ao_basis_ao_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*ao_num*5_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in get_ao_vgl' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + integer :: i,k + do k=1,n_points_final_grid + do i=1,ao_num + aos_in_r_array_qmckl(i,k) = vgl(i,1,k) + aos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k) + aos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k) + aos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k) + aos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, mos_in_r_array_qmckl, (mo_num,n_points_final_grid)] +&BEGIN_PROVIDER [ double precision, mos_grad_in_r_array_qmckl, (mo_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [ double precision, mos_lapl_in_r_array_qmckl, (mo_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! moS computed with qmckl + END_DOC + use qmckl + + integer*8 :: qmckl_ctx + integer(qmckl_exit_code) :: rc + + qmckl_ctx = qmckl_context_create() + + rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename))) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in read_trexio' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in set_electron_point' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + double precision, allocatable :: vgl(:,:,:) + allocate( vgl(mo_num,5,n_points_final_grid)) + rc = qmckl_get_mo_basis_mo_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in get_mo_vgl' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + integer :: i,k + do k=1,n_points_final_grid + do i=1,mo_num + mos_in_r_array_qmckl(i,k) = vgl(i,1,k) + mos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k) + mos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k) + mos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k) + mos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k) + enddo + enddo + +END_PROVIDER + +