From 33ff78d7b3649751d88c5d091c6aa769431f8c99 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Jan 2025 16:04:20 +0100 Subject: [PATCH 1/4] Fixed TREXIO norm problems --- ocaml/Input_ao_basis.ml | 44 ------------------------- scripts/qp_import_trexio.py | 10 +++--- src/ao_basis/EZFIO.cfg | 12 ------- src/ao_basis/aos.irp.f | 3 +- src/basis/EZFIO.cfg | 12 +++++++ src/basis/basis.irp.f | 17 ++++++++-- src/trexio/export_trexio_routines.irp.f | 27 ++++++--------- 7 files changed, 44 insertions(+), 81 deletions(-) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index d9e28e048..343a4ae0b 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -13,8 +13,6 @@ module Ao_basis : sig ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - ao_normalized : bool; - primitives_normalized : bool; } [@@deriving sexp] ;; val read : unit -> t option @@ -36,8 +34,6 @@ end = struct ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - ao_normalized : bool; - primitives_normalized : bool; } [@@deriving sexp] ;; @@ -133,24 +129,6 @@ end = struct Ezfio.get_ao_basis_ao_cartesian () ;; - let read_ao_normalized () = - if not (Ezfio.has_ao_basis_ao_normalized()) then - get_default "ao_normalized" - |> bool_of_string - |> Ezfio.set_ao_basis_ao_normalized - ; - Ezfio.get_ao_basis_ao_normalized () - ;; - - let read_primitives_normalized () = - if not (Ezfio.has_ao_basis_primitives_normalized()) then - get_default "primitives_normalized" - |> bool_of_string - |> Ezfio.set_ao_basis_primitives_normalized - ; - Ezfio.get_ao_basis_primitives_normalized () - ;; - let to_long_basis b = let ao_num = AO_number.to_int b.ao_num in let gto_array = Array.init (AO_number.to_int b.ao_num) @@ -213,8 +191,6 @@ end = struct ao_coef ; ao_expo ; ao_cartesian ; - ao_normalized ; - primitives_normalized ; } = b in write_md5 b ; @@ -247,8 +223,6 @@ end = struct ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; Ezfio.set_ao_basis_ao_cartesian(ao_cartesian); - Ezfio.set_ao_basis_ao_normalized(ao_normalized); - Ezfio.set_ao_basis_primitives_normalized(primitives_normalized); let ao_coef = Array.to_list ao_coef @@ -280,8 +254,6 @@ end = struct ao_coef = read_ao_coef () ; ao_expo = read_ao_expo () ; ao_cartesian = read_ao_cartesian () ; - ao_normalized = read_ao_normalized () ; - primitives_normalized = read_primitives_normalized () ; } in to_md5 result @@ -392,8 +364,6 @@ end = struct { ao_basis = name ; ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ; ao_power ; ao_coef ; ao_expo ; ao_cartesian ; - ao_normalized = bool_of_string @@ get_default "ao_normalized"; - primitives_normalized = bool_of_string @@ get_default "primitives_normalized"; } ;; @@ -448,14 +418,6 @@ Cartesian coordinates (6d,10f,...) :: ao_cartesian = %s -Use normalized primitive functions :: - - primitives_normalized = %s - -Use normalized basis functions :: - - ao_normalized = %s - Basis set (read-only) :: %s @@ -469,8 +431,6 @@ Basis set (read-only) :: " (AO_basis_name.to_string b.ao_basis) (string_of_bool b.ao_cartesian) - (string_of_bool b.primitives_normalized) - (string_of_bool b.ao_normalized) (Basis.to_string short_basis |> String_ext.split ~on:'\n' |> list_map (fun x-> " "^x) @@ -507,8 +467,6 @@ ao_power = %s ao_coef = %s ao_expo = %s ao_cartesian = %s -ao_normalized = %s -primitives_normalized = %s md5 = %s " (AO_basis_name.to_string b.ao_basis) @@ -525,8 +483,6 @@ md5 = %s (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string |> String.concat ", ") (b.ao_cartesian |> string_of_bool) - (b.ao_normalized |> string_of_bool) - (b.primitives_normalized |> string_of_bool) (to_md5 b |> MD5.to_string ) ;; diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index fc76f8de3..ef4ac12ac 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -193,9 +193,11 @@ def write_ezfio(trexio_filename, filename): shell_factor = trexio.read_basis_shell_factor(trexio_file) prim_factor = trexio.read_basis_prim_factor(trexio_file) - for i,p in enumerate(prim_factor): - coefficient[i] *= prim_factor[i] - ezfio.set_ao_basis_primitives_normalized(False) + ezfio.set_basis_prim_normalization_factor(prim_factor) + ezfio.set_basis_primitives_normalized(True) + ezfio.set_basis_ao_normalized(False) + for i, shell_idx in enumerate(shell_index): + coefficient[i] *= shell_factor[shell_idx] ezfio.set_basis_prim_coef(coefficient) elif basis_type.lower() == "numerical": @@ -391,7 +393,7 @@ def write_ezfio(trexio_filename, filename): # Renormalize MO coefs if needed if trexio.has_ao_normalization(trexio_file_cart): - ezfio.set_ao_basis_ao_normalized(False) + ezfio.set_basis_ao_normalized(False) norm = trexio.read_ao_normalization(trexio_file_cart) # for j in range(mo_num): # for i,f in enumerate(norm): diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index bd716383c..c22f8029e 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -55,18 +55,6 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...) interface: ezfio, provider default: false -[ao_normalized] -type: logical -doc: Use normalized basis functions -interface: ezfio, provider -default: true - -[primitives_normalized] -type: logical -doc: Use normalized primitive functions -interface: ezfio, provider -default: true - [use_cgtos] type: logical doc: If true, use cgtos for AO integrals diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 348533980..d718e9353 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -82,10 +82,11 @@ enddo ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm) - if (.not.ao_normalized) then + if (ao_normalized) then do j=1,ao_prim_num(i) ao_coef_normalized(i,j) = ao_coef_normalized(i,j) * ao_coef_normalization_factor(i) enddo + else ao_coef_normalization_factor(i) = 1.d0 endif enddo diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index a6864418a..03e224e4c 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -73,3 +73,15 @@ size: (basis.prim_num) interface: ezfio, provider +[primitives_normalized] +type: logical +doc: If true, assume primitive basis functions are normalized +interface: ezfio, provider, ocaml +default: true + +[ao_normalized] +type: logical +doc: If true, normalize the basis functions +interface: ezfio, provider, ocaml +default: false + diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index b750d75a4..5374e5be9 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -6,6 +6,11 @@ logical :: has PROVIDE ezfio_filename + if (.not.ao_normalized) then + shell_normalization_factor = 1.d0 + return + endif + if (mpi_master) then if (size(shell_normalization_factor) == 0) return @@ -70,6 +75,12 @@ logical :: has PROVIDE ezfio_filename + + if (.not.primitives_normalized) then + prim_normalization_factor(:) = 1.d0 + return + endif + if (mpi_master) then if (size(prim_normalization_factor) == 0) return @@ -95,9 +106,9 @@ do k=1, prim_num if (shell_index(k) /= i) cycle - call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & - powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) - prim_normalization_factor(k) = 1.d0/dsqrt(norm) + call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + prim_normalization_factor(k) = 1.d0/dsqrt(norm) enddo enddo diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f index 0774bcd91..c60b1aa00 100644 --- a/src/trexio/export_trexio_routines.irp.f +++ b/src/trexio/export_trexio_routines.irp.f @@ -271,11 +271,7 @@ subroutine export_trexio(update,full_path) call trexio_assert(rc, TREXIO_SUCCESS) allocate(factor(shell_num)) -! if (ao_normalized) then - factor(1:shell_num) = shell_normalization_factor(1:shell_num) -! else -! factor(1:shell_num) = 1.d0 -! endif + factor(1:shell_num) = shell_normalization_factor(1:shell_num) rc = trexio_write_basis_shell_factor(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) @@ -291,11 +287,12 @@ subroutine export_trexio(update,full_path) call trexio_assert(rc, TREXIO_SUCCESS) allocate(factor(prim_num)) -! if (primitives_normalized) then + if (primitives_normalized) then factor(1:prim_num) = prim_normalization_factor(1:prim_num) -! else -! factor(1:prim_num) = 1.d0 -! endif + else + factor(1:prim_num) = 1.d0 + endif + rc = trexio_write_basis_prim_factor(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) deallocate(factor) @@ -324,14 +321,10 @@ subroutine export_trexio(update,full_path) C_A(3) = 0.d0 allocate(factor(ao_num)) - if (ao_normalized) then - do i=1,ao_num - l = ao_first_of_shell(ao_shell(i)) - factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) - enddo - else - factor(:) = 1.d0 - endif + do i=1,ao_num + l = ao_first_of_shell(ao_shell(i)) + factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) + enddo rc = trexio_write_ao_normalization(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) deallocate(factor) From 1a5f0b3e38c4f32f345c547f5bbd3a69af08e0f2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Jan 2025 17:48:37 +0100 Subject: [PATCH 2/4] Accelerated cholesky integrals in pt2 --- src/mo_two_e_ints/map_integrals.irp.f | 48 +++++++++++++++++++++------ 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 9f485b799..a1471a1cd 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -143,7 +143,7 @@ double precision function get_two_e_integral(i,j,k,l,map) END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx - integer :: ii + integer :: ii, kk type(map_type), intent(inout) :: map real(integral_kind) :: tmp @@ -177,10 +177,10 @@ double precision function get_two_e_integral(i,j,k,l,map) if (do_mo_cholesky) then - double precision, external :: ddot - get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1) -! double precision, external :: get_from_mo_cholesky_cache -! get_two_e_integral = get_from_mo_cholesky_cache(i,j,k,l,.False.) + get_two_e_integral = 0.d0 + do kk=1,cholesky_mo_num + get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo else @@ -516,11 +516,39 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) if (do_mo_cholesky) then - double precision, external :: ddot - do i=1,sze - out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & - cholesky_mo_transp(1,i,l), 1) - enddo + if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. & + (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then + + integer :: kk + + do i=1,mo_integrals_cache_min-1 + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + do i=mo_integrals_cache_min,mo_integrals_cache_max + out_val(i) = get_two_e_integral_cache(i,i,k,l) + enddo + + do i=mo_integrals_cache_max, sze + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + else + + do i=1,sze + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + endif else From db443e0a5af02d2fca237f58e3a78bf64560a3cd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Jan 2025 18:20:43 +0100 Subject: [PATCH 3/4] Put back ddot for AMD --- src/mo_two_e_ints/map_integrals.irp.f | 42 +++++++++++++++++---------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index a1471a1cd..76f169b41 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -177,10 +177,13 @@ double precision function get_two_e_integral(i,j,k,l,map) if (do_mo_cholesky) then - get_two_e_integral = 0.d0 - do kk=1,cholesky_mo_num - get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + double precision, external :: ddot + get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1) + +! get_two_e_integral = 0.d0 +! do kk=1,cholesky_mo_num +! get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo else @@ -519,13 +522,16 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. & (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then + double precision, external :: ddot integer :: kk do i=1,mo_integrals_cache_min-1 - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo do i=mo_integrals_cache_min,mo_integrals_cache_max @@ -533,19 +539,23 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) enddo do i=mo_integrals_cache_max, sze - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo else do i=1,sze - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo endif From ef234305e701ae50bf600127631fc9a8d84f4792 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 22 Jan 2025 11:06:53 +0100 Subject: [PATCH 4/4] Fixed ao_normalization. CP2K OK --- scripts/qp_import_trexio.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index ef4ac12ac..23f48eef2 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -191,11 +191,12 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_nucleus_shell_num(nucl_shell_num) - shell_factor = trexio.read_basis_shell_factor(trexio_file) prim_factor = trexio.read_basis_prim_factor(trexio_file) ezfio.set_basis_prim_normalization_factor(prim_factor) ezfio.set_basis_primitives_normalized(True) ezfio.set_basis_ao_normalized(False) + + shell_factor = trexio.read_basis_shell_factor(trexio_file) for i, shell_idx in enumerate(shell_index): coefficient[i] *= shell_factor[shell_idx] ezfio.set_basis_prim_coef(coefficient) @@ -249,7 +250,6 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_shell_index([x+1 for x in shell_index]) ezfio.set_basis_nucleus_shell_num(nucl_shell_num) - shell_factor = trexio.read_basis_shell_factor(trexio_file) else: raise TypeError @@ -319,13 +319,18 @@ def write_ezfio(trexio_filename, filename): exponent.append(expo[i]) num_prim.append(num_prim0[i]) - print (len(coefficient), ao_num) assert (len(coefficient) == ao_num) + ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) ezfio.set_ao_basis_ao_prim_num(num_prim) prim_num_max = max( [ len(x) for x in coefficient ] ) + ao_normalization = trexio.read_ao_normalization(trexio_file_cart) + for i, coef in enumerate(coefficient): + for j in range(len(coef)): + coef[j] *= ao_normalization[i] + for i in range(ao_num): coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] exponent [i] += [0. for j in range(len(exponent[i]), prim_num_max)] @@ -340,7 +345,6 @@ def write_ezfio(trexio_filename, filename): coef.append(coefficient[j]) expo.append(exponent[j]) -# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max) ezfio.set_ao_basis_ao_coef(coef) ezfio.set_ao_basis_ao_expo(expo) @@ -390,14 +394,6 @@ def write_ezfio(trexio_filename, filename): # Read coefs from temporary cartesian file created in the AO section MoMatrix = trexio.read_mo_coefficient(trexio_file_cart) - - # Renormalize MO coefs if needed - if trexio.has_ao_normalization(trexio_file_cart): - ezfio.set_basis_ao_normalized(False) - norm = trexio.read_ao_normalization(trexio_file_cart) -# for j in range(mo_num): -# for i,f in enumerate(norm): -# MoMatrix[i,j] *= f ezfio.set_mo_basis_mo_coef(MoMatrix) mo_occ = [ 0. for i in range(mo_num) ]