Skip to content

Commit

Permalink
Merge pull request #363 from QuantumPackage/fix-trexio
Browse files Browse the repository at this point in the history
Fixed TREXIO norm problems
  • Loading branch information
scemama authored Jan 22, 2025
2 parents 360aae7 + ef23430 commit b9d178e
Showing 8 changed files with 97 additions and 100 deletions.
44 changes: 0 additions & 44 deletions ocaml/Input_ao_basis.ml
Original file line number Diff line number Diff line change
@@ -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 )

;;
28 changes: 13 additions & 15 deletions scripts/qp_import_trexio.py
Original file line number Diff line number Diff line change
@@ -191,11 +191,14 @@ 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)
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)

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)

elif basis_type.lower() == "numerical":
@@ -247,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

@@ -317,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)]
@@ -338,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)

@@ -388,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_ao_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) ]
12 changes: 0 additions & 12 deletions src/ao_basis/EZFIO.cfg
Original file line number Diff line number Diff line change
@@ -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
3 changes: 2 additions & 1 deletion src/ao_basis/aos.irp.f
Original file line number Diff line number Diff line change
@@ -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
12 changes: 12 additions & 0 deletions src/basis/EZFIO.cfg
Original file line number Diff line number Diff line change
@@ -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

17 changes: 14 additions & 3 deletions src/basis/basis.irp.f
Original file line number Diff line number Diff line change
@@ -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

54 changes: 46 additions & 8 deletions src/mo_two_e_ints/map_integrals.irp.f
Original file line number Diff line number Diff line change
@@ -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

@@ -179,8 +179,11 @@ double precision function get_two_e_integral(i,j,k,l,map)

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 +519,46 @@ 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

double precision, external :: ddot
integer :: kk

do i=1,mo_integrals_cache_min-1
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
out_val(i) = get_two_e_integral_cache(i,i,k,l)
enddo

do i=mo_integrals_cache_max, sze
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) = 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

else

Loading

0 comments on commit b9d178e

Please sign in to comment.