From 2cffbdcc9d253f4feb1e0eeb5b13ba1afba77fca Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 12 Feb 2020 16:34:32 -0600 Subject: [PATCH] significant restructuring of complex int parts instead of real/imag parts read separately, use ezfio to read/write complex arrays with extra dimension of size 2 converter needs to be tested (might need to transpose some axes in arrays) converter has extra garbage that needs to be removed after testing --- REPLACE | 3 + ocaml/Input_mo_basis.ml | 86 ++---- ocaml/qptypes_generator.ml | 1 + src/ao_basis/EZFIO.cfg | 2 +- src/ao_basis/aos_complex.irp.f | 4 +- src/ao_one_e_ints/EZFIO.cfg | 30 +- src/ao_one_e_ints/ao_one_e_ints.irp.f | 77 +++-- src/ao_one_e_ints/ao_overlap.irp.f | 48 +-- src/ao_one_e_ints/kin_ao_ints.irp.f | 60 ++-- src/ao_one_e_ints/pot_ao_ints.irp.f | 55 ++-- src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f | 49 +-- src/ao_two_e_ints/EZFIO.cfg | 10 +- src/ao_two_e_ints/df_ao_ints.irp.f | 281 ++++++++++++++++-- src/ao_two_e_ints/map_integrals_complex.irp.f | 20 ++ src/ao_two_e_ints/two_e_integrals.irp.f | 5 + src/mo_basis/EZFIO.cfg | 14 +- src/mo_basis/mos_complex.irp.f | 90 ++++-- src/mo_basis/utils.irp.f | 48 ++- src/mo_two_e_ints/EZFIO.cfg | 4 +- src/mo_two_e_ints/df_mo_ints.irp.f | 16 +- src/scf_utils/huckel_complex.irp.f | 4 +- .../Gen_Ezfio_from_integral_complex_3idx.sh | 6 +- src/utils_complex/MolPyscfToQPkpts.py | 2 +- .../create_ezfio_complex_3idx.py | 70 +++-- 24 files changed, 648 insertions(+), 337 deletions(-) diff --git a/REPLACE b/REPLACE index 42d530b0..8ec6d09d 100755 --- a/REPLACE +++ b/REPLACE @@ -869,3 +869,6 @@ qp_name mo_two_e_integral_periodic -r mo_two_e_integral_complex qp_name get_mo_two_e_integrals_ij_periodic -r get_mo_two_e_integrals_ij_complex qp_name get_mo_two_e_integrals_coulomb_ii_periodic -r get_mo_two_e_integrals_coulomb_ii_complex qp_name get_mo_two_e_integrals_coulomb_ijij_periodic -r get_mo_two_e_integrals_coulomb_ijij_complex +qp_name ao_kpt_num -r ao_num_per_kpt +qp_name mo_kpt_num -r mo_num_per_kpt +qp_name num_kpts -r kpt_num diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index 94435349..46f8240e 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -2,7 +2,6 @@ open Qptypes open Qputils open Sexplib.Std - module Mo_basis : sig type t = { mo_num : MO_number.t ; @@ -10,7 +9,6 @@ module Mo_basis : sig mo_class : MO_class.t array; mo_occ : MO_occ.t array; mo_coef : (MO_coef.t array) array; - mo_coef_imag : (MO_coef.t array) array option; ao_md5 : MD5.t; } [@@deriving sexp] val read : unit -> t option @@ -25,11 +23,13 @@ end = struct mo_class : MO_class.t array; mo_occ : MO_occ.t array; mo_coef : (MO_coef.t array) array; - mo_coef_imag : (MO_coef.t array) array option; ao_md5 : MD5.t; } [@@deriving sexp] + let get_default = Qpackage.get_ezfio_default "mo_basis" + let is_complex = lazy (Ezfio.get_nuclei_is_complex () ) + let read_mo_label () = if not (Ezfio.has_mo_basis_mo_label ()) then Ezfio.set_mo_basis_mo_label "None" @@ -43,14 +43,7 @@ end = struct mo_coef = Array.map (fun mo -> Array.init (Array.length mo) (fun i -> mo.(ordering.(i))) - ) b.mo_coef ; - mo_coef_imag = - match b.mo_coef_imag with - | None -> None - | Some x -> Some ( Array.map (fun mo -> - Array.init (Array.length mo) - (fun i -> mo.(ordering.(i))) - ) x ) + ) b.mo_coef } let read_ao_md5 () = @@ -69,7 +62,10 @@ end = struct |> MD5.of_string in if (ao_md5 <> result) then - failwith "The current MOs don't correspond to the current AOs."; + begin + Printf.eprintf ":%s:\n:%s:\n%!" (MD5.to_string ao_md5) (MD5.to_string result); + failwith "The current MOs don't correspond to the current AOs." + end; result @@ -120,29 +116,21 @@ end = struct let read_mo_coef () = - let a = Ezfio.get_mo_basis_mo_coef () - |> Ezfio.flattened_ezfio - |> Array.map MO_coef.of_float + let a = + ( + if Lazy.force is_complex then + Ezfio.get_mo_basis_mo_coef_complex () + else + Ezfio.get_mo_basis_mo_coef () + ) + |> Ezfio.flattened_ezfio + |> Array.map MO_coef.of_float in let mo_num = read_mo_num () |> MO_number.to_int in let ao_num = (Array.length a)/mo_num in - Array.init mo_num (fun j -> - Array.sub a (j*ao_num) (ao_num) - ) - - let read_mo_coef_imag () = - if Ezfio.has_mo_basis_mo_coef_imag () then - let a = - Ezfio.get_mo_basis_mo_coef_imag () - |> Ezfio.flattened_ezfio - |> Array.map MO_coef.of_float - in - let mo_num = read_mo_num () |> MO_number.to_int in - let ao_num = (Array.length a)/mo_num in - Some (Array.init mo_num (fun j -> + Array.init mo_num (fun j -> Array.sub a (j*ao_num) (ao_num) - ) ) - else None + ) let read () = @@ -153,7 +141,6 @@ end = struct mo_class = read_mo_class (); mo_occ = read_mo_occ (); mo_coef = read_mo_coef (); - mo_coef_imag = read_mo_coef_imag (); ao_md5 = read_ao_md5 (); } else @@ -161,7 +148,6 @@ end = struct let mo_coef_to_string mo_coef = - (*TODO : add imaginary part here *) let ao_num = Array.length mo_coef.(0) and mo_num = Array.length mo_coef in let rec print_five imin imax = @@ -247,7 +233,6 @@ MO coefficients :: let to_string b = - (*TODO : add imaginary part here *) Printf.sprintf " mo_label = \"%s\" mo_num = %s @@ -300,31 +285,22 @@ mo_coef = %s let write_mo_coef a = let mo_num = Array.length a in - let ao_num = Array.length a.(0) in + let ao_num = + let x = Array.length a.(0) in + if Lazy.force is_complex then x/2 else x + in let data = Array.map (fun mo -> Array.map MO_coef.to_float mo |> Array.to_list) a |> Array.to_list |> List.concat - in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data - |> Ezfio.set_mo_basis_mo_coef - - - let write_mo_coef_imag a = - match a with - | None -> () - | Some a -> - begin - let mo_num = Array.length a in - let ao_num = Array.length a.(0) in - let data = - Array.map (fun mo -> Array.map MO_coef.to_float mo - |> Array.to_list) a - |> Array.to_list - |> List.concat - in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data - |> Ezfio.set_mo_basis_mo_coef_imag - end + in + if Lazy.force is_complex then + (Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| 2 ; ao_num ; mo_num |] ~data + |> Ezfio.set_mo_basis_mo_coef_complex ) + else + (Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data + |> Ezfio.set_mo_basis_mo_coef ) let write @@ -333,7 +309,6 @@ mo_coef = %s mo_class : MO_class.t array; mo_occ : MO_occ.t array; mo_coef : (MO_coef.t array) array; - mo_coef_imag : (MO_coef.t array) array option; ao_md5 : MD5.t; } = write_mo_num mo_num; @@ -341,7 +316,6 @@ mo_coef = %s write_mo_class mo_class; write_mo_occ mo_occ; write_mo_coef mo_coef; - write_mo_coef_imag mo_coef_imag; write_md5 ao_md5 diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 2c54a218..610b67d1 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -166,6 +166,7 @@ let input_ezfio = " let untouched = " + module MO_guess : sig type t [@@deriving sexp] val to_string : t -> string diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index 4b17acad..337b1fee 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -55,7 +55,7 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...) interface: ezfio, provider default: false -[ao_kpt_num] +[ao_num_per_kpt] type: integer doc: Number of |AOs| per kpt default: =(ao_basis.ao_num/nuclei.kpt_num) diff --git a/src/ao_basis/aos_complex.irp.f b/src/ao_basis/aos_complex.irp.f index 8ed10c43..afec0548 100644 --- a/src/ao_basis/aos_complex.irp.f +++ b/src/ao_basis/aos_complex.irp.f @@ -1,7 +1,7 @@ -BEGIN_PROVIDER [ integer, ao_kpt_num ] +BEGIN_PROVIDER [ integer, ao_num_per_kpt ] implicit none BEGIN_DOC ! number of aos per kpt. END_DOC - ao_kpt_num = ao_num/kpt_num + ao_num_per_kpt = ao_num/kpt_num END_PROVIDER diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 9ef019fa..583c7757 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -4,10 +4,10 @@ doc: Nucleus-electron integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[ao_integrals_n_e_imag] +[ao_integrals_n_e_complex] type: double precision -doc: Imaginary part of the nucleus-electron integrals in |AO| basis set -size: (ao_basis.ao_num,ao_basis.ao_num) +doc: Complex nucleus-electron integrals in |AO| basis set +size: (2,ao_basis.ao_num,ao_basis.ao_num) interface: ezfio [io_ao_integrals_n_e] @@ -23,10 +23,10 @@ doc: Kinetic energy integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[ao_integrals_kinetic_imag] +[ao_integrals_kinetic_complex] type: double precision -doc: Imaginary part of the kinetic energy integrals in |AO| basis set -size: (ao_basis.ao_num,ao_basis.ao_num) +doc: Complex kinetic energy integrals in |AO| basis set +size: (2,ao_basis.ao_num,ao_basis.ao_num) interface: ezfio [io_ao_integrals_kinetic] @@ -42,10 +42,10 @@ doc: Pseudopotential integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[ao_integrals_pseudo_imag] +[ao_integrals_pseudo_complex] type: double precision -doc: Imaginary part of the pseudopotential integrals in |AO| basis set -size: (ao_basis.ao_num,ao_basis.ao_num) +doc: Complex pseudopotential integrals in |AO| basis set +size: (2,ao_basis.ao_num,ao_basis.ao_num) interface: ezfio [io_ao_integrals_pseudo] @@ -61,10 +61,10 @@ doc: Overlap integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[ao_integrals_overlap_imag] +[ao_integrals_overlap_complex] type: double precision -doc: Imaginary part of the overlap integrals in |AO| basis set -size: (ao_basis.ao_num,ao_basis.ao_num) +doc: Complex overlap integrals in |AO| basis set +size: (2,ao_basis.ao_num,ao_basis.ao_num) interface: ezfio [io_ao_integrals_overlap] @@ -80,10 +80,10 @@ doc: Combined integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[ao_one_e_integrals_imag] +[ao_one_e_integrals_complex] type: double precision -doc: Imaginary part of the combined integrals in |AO| basis set -size: (ao_basis.ao_num,ao_basis.ao_num) +doc: Complex combined integrals in |AO| basis set +size: (2,ao_basis.ao_num,ao_basis.ao_num) interface: ezfio [io_ao_one_e_integrals] diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f index b5e8872e..be70bf23 100644 --- a/src/ao_one_e_ints/ao_one_e_ints.irp.f +++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f @@ -5,7 +5,10 @@ BEGIN_DOC ! One-electron Hamiltonian in the |AO| basis. END_DOC - + if (is_complex) then + print*,"you shouldn't be here for complex",irp_here + stop -1 + endif IF (read_ao_one_e_integrals) THEN call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals) ELSE @@ -27,43 +30,55 @@ END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] - implicit none - integer :: i,j,n,l - BEGIN_DOC - ! One-electron Hamiltonian in the |AO| basis. - END_DOC +!BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] +! implicit none +! integer :: i,j,n,l +! BEGIN_DOC +! ! One-electron Hamiltonian in the |AO| basis. +! END_DOC +! +! IF (read_ao_one_e_integrals) THEN +! call ezfio_get_ao_one_e_ints_ao_one_e_integrals_imag(ao_one_e_integrals_imag) +! ELSE +! ao_one_e_integrals_imag = ao_integrals_n_e_imag + ao_kinetic_integrals_imag +! +! IF (DO_PSEUDO) THEN +! ao_one_e_integrals_imag += ao_pseudo_integrals_imag +! ENDIF +! ENDIF +! +! IF (write_ao_one_e_integrals) THEN +! call ezfio_set_ao_one_e_ints_ao_one_e_integrals_imag(ao_one_e_integrals_imag) +! print *, 'AO one-e integrals written to disk' +! ENDIF +! +!END_PROVIDER - IF (read_ao_one_e_integrals) THEN - call ezfio_get_ao_one_e_ints_ao_one_e_integrals_imag(ao_one_e_integrals_imag) - ELSE - ao_one_e_integrals_imag = ao_integrals_n_e_imag + ao_kinetic_integrals_imag - - IF (DO_PSEUDO) THEN - ao_one_e_integrals_imag += ao_pseudo_integrals_imag - ENDIF - ENDIF - - IF (write_ao_one_e_integrals) THEN - call ezfio_set_ao_one_e_ints_ao_one_e_integrals_imag(ao_one_e_integrals_imag) - print *, 'AO one-e integrals written to disk' - ENDIF - -END_PROVIDER - -BEGIN_PROVIDER [ complex*16, ao_one_e_integrals_complex,(ao_num,ao_num)] + BEGIN_PROVIDER [ complex*16, ao_one_e_integrals_complex,(ao_num,ao_num)] +&BEGIN_PROVIDER [ double precision, ao_one_e_integrals_diag_complex,(ao_num)] implicit none integer :: i,j,n,l BEGIN_DOC ! One-electron Hamiltonian in the |AO| basis. END_DOC - do i=1,ao_num - do j=1,ao_num - ao_one_e_integrals_complex(j,i)=dcmplx(ao_one_e_integrals(j,i), & - ao_one_e_integrals_imag(j,i)) - enddo - enddo + IF (read_ao_one_e_integrals) THEN + call ezfio_get_ao_one_e_ints_ao_one_e_integrals_complex(ao_one_e_integrals_complex) + ELSE + ao_one_e_integrals_complex = ao_integrals_n_e_complex + ao_kinetic_integrals_complex + IF (DO_PSEUDO) THEN + ao_one_e_integrals_complex += ao_pseudo_integrals_complex + ENDIF + ENDIF + + DO j = 1, ao_num + ao_one_e_integrals_diag_complex(j) = dble(ao_one_e_integrals_complex(j,j)) + ENDDO + + IF (write_ao_one_e_integrals) THEN + call ezfio_set_ao_one_e_ints_ao_one_e_integrals_complex(ao_one_e_integrals_complex) + print *, 'AO one-e integrals written to disk' + ENDIF END_PROVIDER diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 52a0ea1c..5afadbbe 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -70,34 +70,38 @@ END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Imaginary part of the overlap - END_DOC - if (read_ao_integrals_overlap) then - call ezfio_get_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) - print *, 'AO overlap integrals read from disk' - else - ao_overlap_imag = 0.d0 - endif - if (write_ao_integrals_overlap) then - call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) - print *, 'AO overlap integrals written to disk' - endif -END_PROVIDER +!BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ] +! implicit none +! BEGIN_DOC +! ! Imaginary part of the overlap +! END_DOC +! if (read_ao_integrals_overlap) then +! call ezfio_get_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) +! print *, 'AO overlap integrals read from disk' +! else +! ao_overlap_imag = 0.d0 +! endif +! if (write_ao_integrals_overlap) then +! call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) +! print *, 'AO overlap integrals written to disk' +! endif +!END_PROVIDER BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ] implicit none BEGIN_DOC ! Overlap for complex AOs END_DOC - integer :: i,j - do j=1,ao_num - do i=1,ao_num - ao_overlap_complex(i,j) = dcmplx( ao_overlap(i,j), ao_overlap_imag(i,j) ) - enddo - enddo + if (read_ao_integrals_overlap) then + call ezfio_get_ao_one_e_ints_ao_integrals_overlap_complex(ao_overlap_complex) + print *, 'AO overlap integrals read from disk' + else + print*,'complex AO overlap ints must be provided',irp_here + endif + if (write_ao_integrals_overlap) then + call ezfio_set_ao_one_e_ints_ao_integrals_overlap_complex(ao_overlap_complex) + print *, 'AO overlap integrals written to disk' + endif END_PROVIDER diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f index ca50114c..f352d1c4 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -149,27 +149,27 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)] endif END_PROVIDER -BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)] - implicit none - BEGIN_DOC - ! Kinetic energy integrals in the |AO| basis. - ! - ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ - ! - END_DOC - integer :: i,j,k,l - - if (read_ao_integrals_kinetic) then - call ezfio_get_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) - print *, 'AO kinetic integrals read from disk' - else - print *, irp_here, ': Not yet implemented' - endif - if (write_ao_integrals_kinetic) then - call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) - print *, 'AO kinetic integrals written to disk' - endif -END_PROVIDER +!BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)] +! implicit none +! BEGIN_DOC +! ! Kinetic energy integrals in the |AO| basis. +! ! +! ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ +! ! +! END_DOC +! integer :: i,j,k,l +! +! if (read_ao_integrals_kinetic) then +! call ezfio_get_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) +! print *, 'AO kinetic integrals read from disk' +! else +! print *, irp_here, ': Not yet implemented' +! endif +! if (write_ao_integrals_kinetic) then +! call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) +! print *, 'AO kinetic integrals written to disk' +! endif +!END_PROVIDER BEGIN_PROVIDER [complex*16, ao_kinetic_integrals_complex, (ao_num,ao_num)] implicit none @@ -179,11 +179,15 @@ BEGIN_PROVIDER [complex*16, ao_kinetic_integrals_complex, (ao_num,ao_num)] ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ ! END_DOC - integer :: i,j - do i=1,ao_num - do j=1,ao_num - ao_kinetic_integrals_complex(j,i) = dcmplx(ao_kinetic_integrals(j,i), & - ao_kinetic_integrals_imag(j,i)) - enddo - enddo + if (read_ao_integrals_kinetic) then + call ezfio_get_ao_one_e_ints_ao_integrals_kinetic_complex(ao_kinetic_integrals_complex) + print *, 'AO kinetic integrals read from disk' + else + print *, irp_here, ': Not yet implemented' + stop -1 + endif + if (write_ao_integrals_kinetic) then + call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_complex(ao_kinetic_integrals_complex) + print *, 'AO kinetic integrals written to disk' + endif END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 63c02dd2..08c78464 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -83,27 +83,27 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)] - BEGIN_DOC - ! Nucleus-electron interaction, in the |AO| basis set. - ! - ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` - END_DOC - implicit none - double precision :: alpha, beta, gama, delta - integer :: num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m - double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - - if (read_ao_integrals_n_e) then - call ezfio_get_ao_one_e_ints_ao_integrals_n_e_imag(ao_integrals_n_e_imag) - print *, 'AO N-e integrals read from disk' - else - print *, irp_here, ': Not yet implemented' - endif -END_PROVIDER +!BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)] +! BEGIN_DOC +! ! Nucleus-electron interaction, in the |AO| basis set. +! ! +! ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` +! END_DOC +! implicit none +! double precision :: alpha, beta, gama, delta +! integer :: num_A,num_B +! double precision :: A_center(3),B_center(3),C_center(3) +! integer :: power_A(3),power_B(3) +! integer :: i,j,k,l,n_pt_in,m +! double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult +! +! if (read_ao_integrals_n_e) then +! call ezfio_get_ao_one_e_ints_ao_integrals_n_e_imag(ao_integrals_n_e_imag) +! print *, 'AO N-e integrals read from disk' +! else +! print *, irp_here, ': Not yet implemented' +! endif +!END_PROVIDER BEGIN_PROVIDER [complex*16, ao_integrals_n_e_complex, (ao_num,ao_num)] implicit none @@ -112,13 +112,12 @@ BEGIN_PROVIDER [complex*16, ao_integrals_n_e_complex, (ao_num,ao_num)] ! ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` END_DOC - integer :: i,j - do i=1,ao_num - do j=1,ao_num - ao_integrals_n_e_complex(j,i) = dcmplx(ao_integrals_n_e(j,i), & - ao_integrals_n_e_imag(j,i)) - enddo - enddo + if (read_ao_integrals_n_e) then + call ezfio_get_ao_one_e_ints_ao_integrals_n_e_complex(ao_integrals_n_e_complex) + print *, 'AO N-e integrals read from disk' + else + print *, irp_here, ': Not yet implemented' + endif END_PROVIDER BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)] diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index a92ba1f4..0032b2ae 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -27,34 +27,39 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)] END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_imag, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Imaginary part of the pseudo_integrals - END_DOC - if (read_ao_integrals_pseudo) then - call ezfio_get_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) - print *, 'AO pseudo_integrals integrals read from disk' - else - ao_pseudo_integrals_imag = 0.d0 - endif - if (write_ao_integrals_pseudo) then - call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) - print *, 'AO pseudo_integrals integrals written to disk' - endif -END_PROVIDER +!BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_imag, (ao_num, ao_num) ] +! implicit none +! BEGIN_DOC +! ! Imaginary part of the pseudo_integrals +! END_DOC +! if (read_ao_integrals_pseudo) then +! call ezfio_get_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) +! print *, 'AO pseudo_integrals integrals read from disk' +! else +! ao_pseudo_integrals_imag = 0.d0 +! endif +! if (write_ao_integrals_pseudo) then +! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) +! print *, 'AO pseudo_integrals integrals written to disk' +! endif +!END_PROVIDER BEGIN_PROVIDER [ complex*16, ao_pseudo_integrals_complex, (ao_num, ao_num) ] implicit none BEGIN_DOC ! Overlap for complex AOs END_DOC - integer :: i,j - do j=1,ao_num - do i=1,ao_num - ao_pseudo_integrals_complex(i,j) = dcmplx( ao_pseudo_integrals(i,j), ao_pseudo_integrals_imag(i,j) ) - enddo - enddo + if (read_ao_integrals_pseudo) then + call ezfio_get_ao_one_e_ints_ao_integrals_pseudo_complex(ao_pseudo_integrals_complex) + print *, 'AO pseudo_integrals integrals read from disk' + else + print*,irp_here,'not implemented' + stop -1 + endif + if (write_ao_integrals_pseudo) then + call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_complex(ao_pseudo_integrals_complex) + print *, 'AO pseudo_integrals integrals written to disk' + endif END_PROVIDER BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index fad4c836..5b50f718 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -29,15 +29,9 @@ doc: Read/Write df |AO| integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[df_ao_integrals_real] +[df_ao_integrals_complex] type: double precision doc: Real part of the df integrals over AOs -size: (ao_basis.ao_kpt_num,ao_basis.ao_kpt_num,ao_two_e_ints.df_num,nuclei.kpt_pair_num) -interface: ezfio - -[df_ao_integrals_imag] -type: double precision -doc: Imaginary part of the df integrals over AOs -size: (ao_basis.ao_kpt_num,ao_basis.ao_kpt_num,ao_two_e_ints.df_num,nuclei.kpt_pair_num) +size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio diff --git a/src/ao_two_e_ints/df_ao_ints.irp.f b/src/ao_two_e_ints/df_ao_ints.irp.f index 20f5f21c..dd936519 100644 --- a/src/ao_two_e_ints/df_ao_ints.irp.f +++ b/src/ao_two_e_ints/df_ao_ints.irp.f @@ -1,6 +1,52 @@ - BEGIN_PROVIDER [double precision, df_ao_integrals_real, (ao_kpt_num,ao_kpt_num,df_num,kpt_pair_num)] -&BEGIN_PROVIDER [double precision, df_ao_integrals_imag, (ao_kpt_num,ao_kpt_num,df_num,kpt_pair_num)] -&BEGIN_PROVIDER [complex*16, df_ao_integrals_complex, (ao_kpt_num,ao_kpt_num,df_num,kpt_pair_num)] +! BEGIN_PROVIDER [double precision, df_ao_integrals_real, (ao_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num)] +! &BEGIN_PROVIDER [double precision, df_ao_integrals_imag, (ao_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num)] +! &BEGIN_PROVIDER [complex*16, df_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num)] +! implicit none +! BEGIN_DOC +! ! df AO integrals +! END_DOC +! integer :: i,j,k,l +! +! if (read_df_ao_integrals) then +! df_ao_integrals_real = 0.d0 +! df_ao_integrals_imag = 0.d0 +! call ezfio_get_ao_two_e_ints_df_ao_integrals_real(df_ao_integrals_real) +! call ezfio_get_ao_two_e_ints_df_ao_integrals_imag(df_ao_integrals_imag) +! print *, 'df AO integrals read from disk' +! do l=1,kpt_pair_num +! do k=1,df_num +! do j=1,ao_num_per_kpt +! do i=1,ao_num_per_kpt +! df_ao_integrals_complex(i,j,k,l) = dcmplx(df_ao_integrals_real(i,j,k,l), & +! df_ao_integrals_imag(i,j,k,l)) +! enddo +! enddo +! enddo +! enddo +! else +! print*,'df ao integrals must be provided',irp_here +! stop -1 +! endif +! +! if (write_df_ao_integrals) then +! ! this probably shouldn't happen +! do l=1,kpt_pair_num +! do k=1,df_num +! do j=1,ao_num_per_kpt +! do i=1,ao_num_per_kpt +! df_ao_integrals_real(i,j,k,l) = dble(df_ao_integrals_complex(i,j,k,l)) +! df_ao_integrals_imag(i,j,k,l) = dimag(df_ao_integrals_complex(i,j,k,l)) +! enddo +! enddo +! enddo +! enddo +! call ezfio_set_ao_two_e_ints_df_ao_integrals_real(df_ao_integrals_real) +! call ezfio_set_ao_two_e_ints_df_ao_integrals_imag(df_ao_integrals_imag) +! print *, 'df AO integrals written to disk' +! endif +! +! END_PROVIDER + BEGIN_PROVIDER [complex*16, df_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num)] implicit none BEGIN_DOC ! df AO integrals @@ -8,42 +54,217 @@ integer :: i,j,k,l if (read_df_ao_integrals) then - df_ao_integrals_real = 0.d0 - df_ao_integrals_imag = 0.d0 - call ezfio_get_ao_two_e_ints_df_ao_integrals_real(df_ao_integrals_real) - call ezfio_get_ao_two_e_ints_df_ao_integrals_imag(df_ao_integrals_imag) + call ezfio_get_ao_two_e_ints_df_ao_integrals_complex(df_ao_integrals_complex) print *, 'df AO integrals read from disk' - do l=1,kpt_pair_num - do k=1,df_num - do j=1,ao_kpt_num - do i=1,ao_kpt_num - df_ao_integrals_complex(i,j,k,l) = dcmplx(df_ao_integrals_real(i,j,k,l), & - df_ao_integrals_imag(i,j,k,l)) - enddo - enddo - enddo - enddo else print*,'df ao integrals must be provided',irp_here stop -1 endif if (write_df_ao_integrals) then - ! this probably shouldn't happen - do l=1,kpt_pair_num - do k=1,df_num - do j=1,ao_kpt_num - do i=1,ao_kpt_num - df_ao_integrals_real(i,j,k,l) = dble(df_ao_integrals_complex(i,j,k,l)) - df_ao_integrals_imag(i,j,k,l) = dimag(df_ao_integrals_complex(i,j,k,l)) - enddo - enddo - enddo - enddo - call ezfio_set_ao_two_e_ints_df_ao_integrals_real(df_ao_integrals_real) - call ezfio_set_ao_two_e_ints_df_ao_integrals_imag(df_ao_integrals_imag) + call ezfio_set_ao_two_e_ints_df_ao_integrals_complex(df_ao_integrals_complex) print *, 'df AO integrals written to disk' endif END_PROVIDER + +subroutine ao_map_fill_from_df + use map_module + implicit none + BEGIN_DOC + ! fill ao bielec integral map using 3-index df integrals + END_DOC + + integer :: i,k,j,l + integer :: ki,kk,kj,kl + integer :: ii,ik,ij,il + integer :: kikk2,kjkl2,jl2,ik2 + integer :: i_ao,j_ao,i_df + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) + + complex*16 :: integral + integer :: n_integrals_1, n_integrals_2 + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:) + double precision :: tmp_re,tmp_im + integer :: ao_num_kpt_2 + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: map_mb + + logical :: use_map1 + integer(keY_kind) :: idx_tmp + double precision :: sign + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) + print*, 'Providing the ao_bielec integrals from 3-index df integrals' + call write_time(6) +! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write') +! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals + + call wall_time(wall_1) + call cpu_time(cpu_1) + + allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,df_num)) + + wall_0 = wall_1 + do kl=1, kpt_num + do kj=1, kl + call idx2_tri_int(kj,kl,kjkl2) + ints_jl = df_ao_integrals_complex(:,:,:,kjkl2) + + !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, & + !$OMP ints_ik, ints_ikjl, i_ao, j_ao, i_df, & + !$OMP n_integrals_1, buffer_i_1, buffer_values_1, & + !$OMP n_integrals_2, buffer_i_2, buffer_values_2, & + !$OMP idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer, kpt_num, df_num, ao_num_per_kpt, ao_num_kpt_2, & + !$OMP kl,kj,kjkl2,ints_jl, & + !$OMP kconserv, df_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) + + allocate( & + ints_ik(ao_num_per_kpt,ao_num_per_kpt,df_num), & + ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), & + buffer_i_1(size_buffer), & + buffer_i_2(size_buffer), & + buffer_values_1(size_buffer), & + buffer_values_2(size_buffer) & + ) + + !$OMP DO SCHEDULE(guided) + do kk=1,kl + ki=kconserv(kl,kk,kj) + if ((kl == kj) .and. (ki > kk)) cycle + call idx2_tri_int(ki,kk,kikk2) + if (kikk2 > kjkl2) cycle + if (ki >= kk) then + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2)) + enddo + enddo + enddo +! ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/))) + else + ints_ik = df_ao_integrals_complex(:,:,:,kikk2) + endif + + call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, df_num, & + (1.d0,0.d0), ints_ik, ao_num_kpt_2, & + ints_jl, ao_num_kpt_2, & + (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) + + n_integrals_1=0 + n_integrals_2=0 + do il=1,ao_num_per_kpt + l=il+(kl-1)*ao_num_per_kpt + do ij=1,ao_num_per_kpt + j=ij+(kj-1)*ao_num_per_kpt + if (j>l) exit + call idx2_tri_int(j,l,jl2) + do ik=1,ao_num_per_kpt + k=ik+(kk-1)*ao_num_per_kpt + if (k>l) exit + do ii=1,ao_num_per_kpt + i=ii+(ki-1)*ao_num_per_kpt + if ((j==l) .and. (i>k)) exit + call idx2_tri_int(i,k,ik2) + if (ik2 > jl2) exit + integral = ints_ikjl(ii,ik,ij,il) +! print*,i,k,j,l,real(integral),imag(integral) + if (cdabs(integral) < ao_integrals_threshold) then + cycle + endif + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + tmp_re = dble(integral) + tmp_im = dimag(integral) + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re + if (sign.ne.0.d0) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp+1 + buffer_values_1(n_integrals_1)=tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + n_integrals_1 = 0 + endif + else + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp + buffer_values_2(n_integrals_2)=tmp_re + if (sign.ne.0.d0) then + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp+1 + buffer_values_2(n_integrals_2)=tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + n_integrals_2 = 0 + endif + endif + + enddo !ii + enddo !ik + enddo !ij + enddo !il + + if (n_integrals_1 > 0) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + endif + enddo !kk + !$OMP END DO NOWAIT + deallocate( & + ints_ik, & + ints_ikjl, & + buffer_i_1, & + buffer_i_2, & + buffer_values_1, & + buffer_values_2 & + ) + !$OMP END PARALLEL + enddo !kj + call wall_time(wall_2) + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(kl)/float(kpt_num), '% in ', & + wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + endif + + enddo !kl + deallocate( ints_jl ) + + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) + call map_sort(ao_integrals_map_2) + call map_unique(ao_integrals_map_2) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + print*,'AO integrals provided:' + print*,' Size of AO map ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + print*,' Number of AO integrals: ', ao_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end subroutine ao_map_fill_from_df + diff --git a/src/ao_two_e_ints/map_integrals_complex.irp.f b/src/ao_two_e_ints/map_integrals_complex.irp.f index dc4e5542..449bca02 100644 --- a/src/ao_two_e_ints/map_integrals_complex.irp.f +++ b/src/ao_two_e_ints/map_integrals_complex.irp.f @@ -1,6 +1,26 @@ use map_module +subroutine idx2_tri_int(i,j,ij) + implicit none + integer, intent(in) :: i,j + integer, intent(out) :: ij + integer :: p,q + p = max(i,j) + q = min(i,j) + ij = q+ishft(p*p-p,-1) +end + +subroutine idx2_tri_key(i,j,ij) + use map_module + implicit none + integer, intent(in) :: i,j + integer(key_kind), intent(out) :: ij + integer(key_kind) :: p,q + p = max(i,j) + q = min(i,j) + ij = q+ishft(p*p-p,-1) +end subroutine two_e_integrals_index_complex(i,j,k,l,i1,p,q) use map_module implicit none diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index f3bb5c20..8cfd29fb 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -359,6 +359,11 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] print*, 'AO integrals provided (periodic)' ao_two_e_integrals_in_map = .True. return + else if (read_df_ao_integrals) then + call ao_map_fill_from_df + print*, 'AO integrals provided from 3-index ao ints (periodic)' + ao_two_e_integrals_in_map = .True. + return else print*,'calculation of periodic AOs not implemented' stop -1 diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index f667f04f..fd9303aa 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -9,17 +9,11 @@ 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_real] +[mo_coef_complex] type: double precision -doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO| +doc: Complex MO 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| -interface: ezfio -size: (ao_basis.ao_num,mo_basis.mo_num) +size: (2,ao_basis.ao_num,mo_basis.mo_num) [mo_label] type: character*(64) @@ -43,7 +37,7 @@ type: character*(32) doc: MD5 checksum characterizing the |AO| basis set. interface: ezfio -[mo_kpt_num] +[mo_num_per_kpt] type: integer doc: Number of |MOs| per kpt interface: ezfio diff --git a/src/mo_basis/mos_complex.irp.f b/src/mo_basis/mos_complex.irp.f index 35987220..2e2f0786 100644 --- a/src/mo_basis/mos_complex.irp.f +++ b/src/mo_basis/mos_complex.irp.f @@ -1,11 +1,74 @@ -BEGIN_PROVIDER [ integer, mo_kpt_num ] +BEGIN_PROVIDER [ integer, mo_num_per_kpt ] implicit none BEGIN_DOC ! number of mos per kpt. END_DOC - mo_kpt_num = mo_num/kpt_num + mo_num_per_kpt = mo_num/kpt_num END_PROVIDER +!BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] +! implicit none +! BEGIN_DOC +! ! Molecular orbital coefficients on |AO| basis set +! ! +! ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO| +! ! +! ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) +! END_DOC +! integer :: i, j +! double precision, allocatable :: buffer_re(:,:),buffer_im(:,:) +! logical :: exists_re,exists_im,exists +! PROVIDE ezfio_filename +! +! +! if (mpi_master) then +! ! Coefs +! call ezfio_has_mo_basis_mo_coef_real(exists_re) +! call ezfio_has_mo_basis_mo_coef_imag(exists_im) +! exists = (exists_re.and.exists_im) +! 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_real/imag with MPI' +! endif +! IRP_ENDIF +! +! if (exists) then +! if (mpi_master) then +! allocate(buffer_re(ao_num,mo_num),buffer_im(ao_num,mo_num)) +! call ezfio_get_mo_basis_mo_coef_real(buffer_re) +! call ezfio_get_mo_basis_mo_coef_imag(buffer_im) +! write(*,*) 'Read mo_coef_real/imag' +! do i=1,mo_num +! do j=1,ao_num +! mo_coef_complex(j,i) = dcmplx(buffer_re(j,i),buffer_im(j,i)) +! enddo +! enddo +! deallocate(buffer_re,buffer_im) +! endif +! IRP_IF MPI +! call MPI_BCAST( mo_coef_complex, mo_num*ao_num, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD, ierr) +! if (ierr /= MPI_SUCCESS) then +! stop 'Unable to read mo_coef_real with MPI' +! endif +! IRP_ENDIF +! else +! ! Orthonormalized AO basis +! do i=1,mo_num +! do j=1,ao_num +! mo_coef_complex(j,i) = ao_ortho_canonical_coef_complex(j,i) +! enddo +! enddo +! endif +!END_PROVIDER + BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] implicit none BEGIN_DOC @@ -16,16 +79,13 @@ BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) END_DOC integer :: i, j - double precision, allocatable :: buffer_re(:,:),buffer_im(:,:) - logical :: exists_re,exists_im,exists + logical :: exists PROVIDE ezfio_filename if (mpi_master) then ! Coefs - call ezfio_has_mo_basis_mo_coef_real(exists_re) - call ezfio_has_mo_basis_mo_coef_imag(exists_im) - exists = (exists_re.and.exists_im) + call ezfio_has_mo_basis_mo_coef_complex(exists) endif IRP_IF MPI_DEBUG print *, irp_here, mpi_rank @@ -36,27 +96,19 @@ BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] 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_real/imag with MPI' + stop 'Unable to read mo_coef_complex with MPI' endif IRP_ENDIF if (exists) then if (mpi_master) then - allocate(buffer_re(ao_num,mo_num),buffer_im(ao_num,mo_num)) - call ezfio_get_mo_basis_mo_coef_real(buffer_re) - call ezfio_get_mo_basis_mo_coef_imag(buffer_im) - write(*,*) 'Read mo_coef_real/imag' - do i=1,mo_num - do j=1,ao_num - mo_coef_complex(j,i) = dcmplx(buffer_re(j,i),buffer_im(j,i)) - enddo - enddo - deallocate(buffer_re,buffer_im) + call ezfio_get_mo_basis_mo_coef_complex(mo_coef_complex) + write(*,*) 'Read mo_coef_complex' endif IRP_IF MPI call MPI_BCAST( mo_coef_complex, mo_num*ao_num, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then - stop 'Unable to read mo_coef_real with MPI' + stop 'Unable to read mo_coef_complex with MPI' endif IRP_ENDIF else diff --git a/src/mo_basis/utils.irp.f b/src/mo_basis/utils.irp.f index a84f9fb7..5f93bb2f 100644 --- a/src/mo_basis/utils.irp.f +++ b/src/mo_basis/utils.irp.f @@ -1,6 +1,7 @@ subroutine save_mos implicit none - double precision, allocatable :: buffer(:,:),buffer_im(:,:) + double precision, allocatable :: buffer(:,:) + complex*16, allocatable :: buffer_c(:,:) integer :: i,j !TODO: change this for periodic? ! save real/imag parts of mo_coef_complex @@ -11,18 +12,15 @@ subroutine save_mos call ezfio_set_mo_basis_mo_label(mo_label) call ezfio_set_mo_basis_ao_md5(ao_md5) if (is_complex) then - allocate ( buffer(ao_num,mo_num),buffer_im(ao_num,mo_num)) - buffer = 0.d0 - buffer_im = 0.d0 + allocate ( buffer_c(ao_num,mo_num)) + buffer_c = (0.d0,0.d0) do j = 1, mo_num do i = 1, ao_num - buffer(i,j) = dble(mo_coef_complex(i,j)) - buffer_im(i,j) = dimag(mo_coef_complex(i,j)) + buffer_c(i,j) = mo_coef_complex(i,j) enddo enddo - call ezfio_set_mo_basis_mo_coef_real(buffer) - call ezfio_set_mo_basis_mo_coef_imag(buffer_im) - deallocate (buffer,buffer_im) + call ezfio_set_mo_basis_mo_coef_complex(buffer_c) + deallocate (buffer_c) else allocate ( buffer(ao_num,mo_num) ) buffer = 0.d0 @@ -42,7 +40,8 @@ end subroutine save_mos_no_occ implicit none - double precision, allocatable :: buffer(:,:),buffer_im(:,:) + double precision, allocatable :: buffer(:,:) + complex*16, allocatable :: buffer_c(:,:) integer :: i,j call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) @@ -50,18 +49,15 @@ subroutine save_mos_no_occ !call ezfio_set_mo_basis_mo_label(mo_label) !call ezfio_set_mo_basis_ao_md5(ao_md5) if (is_complex) then - allocate ( buffer(ao_num,mo_num),buffer_im(ao_num,mo_num)) - buffer = 0.d0 - buffer_im = 0.d0 + allocate ( buffer_c(ao_num,mo_num)) + buffer_c = (0.d0,0.d0) do j = 1, mo_num do i = 1, ao_num - buffer(i,j) = dble(mo_coef_complex(i,j)) - buffer_im(i,j) = dimag(mo_coef_complex(i,j)) + buffer_c(i,j) = mo_coef_complex(i,j) enddo enddo - call ezfio_set_mo_basis_mo_coef_real(buffer) - call ezfio_set_mo_basis_mo_coef_imag(buffer_im) - deallocate (buffer,buffer_im) + call ezfio_set_mo_basis_mo_coef_complex(buffer_c) + deallocate (buffer_c) else allocate ( buffer(ao_num,mo_num) ) buffer = 0.d0 @@ -78,7 +74,8 @@ end subroutine save_mos_truncated(n) implicit none - double precision, allocatable :: buffer(:,:),buffer_im(:,:) + double precision, allocatable :: buffer(:,:) + complex*16, allocatable :: buffer_c(:,:) integer :: i,j,n call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) @@ -87,18 +84,15 @@ subroutine save_mos_truncated(n) call ezfio_set_mo_basis_mo_label(mo_label) call ezfio_set_mo_basis_ao_md5(ao_md5) if (is_complex) then - allocate ( buffer(ao_num,n),buffer_im(ao_num,n)) - buffer = 0.d0 - buffer_im = 0.d0 + allocate ( buffer_c(ao_num,mo_num)) + buffer_c = (0.d0,0.d0) do j = 1, n do i = 1, ao_num - buffer(i,j) = dble(mo_coef_complex(i,j)) - buffer_im(i,j) = dimag(mo_coef_complex(i,j)) + buffer_c(i,j) = mo_coef_complex(i,j) enddo enddo - call ezfio_set_mo_basis_mo_coef_real(buffer) - call ezfio_set_mo_basis_mo_coef_imag(buffer_im) - deallocate (buffer,buffer_im) + call ezfio_set_mo_basis_mo_coef_complex(buffer_c) + deallocate (buffer_c) else allocate ( buffer(ao_num,n) ) buffer = 0.d0 diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index 92bc086c..fc1ff2e1 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -20,12 +20,12 @@ default: None [df_mo_integrals_real] type: double precision doc: Real part of the df integrals over MOs -size: (mo_basis.mo_kpt_num,mo_basis.mo_kpt_num,ao_two_e_ints.df_num,nuclei.kpt_pair_num) +size: (mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio [df_mo_integrals_imag] type: double precision doc: Imaginary part of the df integrals over MOs -size: (mo_basis.mo_kpt_num,mo_basis.mo_kpt_num,ao_two_e_ints.df_num,nuclei.kpt_pair_num) +size: (mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index 62eb683f..5d85056b 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -1,6 +1,6 @@ - BEGIN_PROVIDER [double precision, df_mo_integrals_real, (mo_kpt_num,mo_kpt_num,df_num,kpt_pair_num)] -&BEGIN_PROVIDER [double precision, df_mo_integrals_imag, (mo_kpt_num,mo_kpt_num,df_num,kpt_pair_num)] -&BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_kpt_num,mo_kpt_num,df_num,kpt_pair_num)] + BEGIN_PROVIDER [double precision, df_mo_integrals_real, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)] +&BEGIN_PROVIDER [double precision, df_mo_integrals_imag, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)] +&BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)] implicit none BEGIN_DOC ! df AO integrals @@ -15,8 +15,8 @@ print *, 'df AO integrals read from disk' do l=1,kpt_pair_num do k=1,df_num - do j=1,mo_kpt_num - do i=1,mo_kpt_num + do j=1,mo_num_per_kpt + do i=1,mo_num_per_kpt df_mo_integrals_complex(i,j,k,l) = dcmplx(df_mo_integrals_real(i,j,k,l), & df_mo_integrals_imag(i,j,k,l)) enddo @@ -24,14 +24,14 @@ enddo enddo else - call df_mo_from_df_ao(df_mo_integrals_complex,df_ao_integrals_complex,mo_kpt_num,ao_kpt_num,df_num,kpt_pair_num) + call df_mo_from_df_ao(df_mo_integrals_complex,df_ao_integrals_complex,mo_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num) endif if (write_df_mo_integrals) then do l=1,kpt_pair_num do k=1,df_num - do j=1,mo_kpt_num - do i=1,mo_kpt_num + do j=1,mo_num_per_kpt + do i=1,mo_num_per_kpt df_mo_integrals_real(i,j,k,l) = dble(df_mo_integrals_complex(i,j,k,l)) df_mo_integrals_imag(i,j,k,l) = dimag(df_mo_integrals_complex(i,j,k,l)) enddo diff --git a/src/scf_utils/huckel_complex.irp.f b/src/scf_utils/huckel_complex.irp.f index 8f448d42..d6da7ffb 100644 --- a/src/scf_utils/huckel_complex.irp.f +++ b/src/scf_utils/huckel_complex.irp.f @@ -15,9 +15,9 @@ subroutine huckel_guess_complex A = 0.d0 do j=1,ao_num do i=1,ao_num - A(i,j) = c * ao_overlap_complex(i,j) * (ao_one_e_integrals_diag(i) + ao_one_e_integrals_diag(j)) + A(i,j) = c * ao_overlap_complex(i,j) * (ao_one_e_integrals_diag_complex(i) + ao_one_e_integrals_diag_complex(j)) enddo - A(j,j) = ao_one_e_integrals_diag(j) + dble(ao_two_e_integral_alpha_complex(j,j)) + A(j,j) = ao_one_e_integrals_diag_complex(j) + dble(ao_two_e_integral_alpha_complex(j,j)) if (dabs(dimag(ao_two_e_integral_alpha_complex(j,j))) .gt. 1.0d-10) then stop 'diagonal elements of ao_bi_elec_integral_alpha should be real' endif diff --git a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh index e560ae38..31d273c1 100755 --- a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh +++ b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh @@ -9,16 +9,16 @@ echo 'Create EZFIO' #read nel nmo natom <<< $(cat param) #read e_nucl <<< $(cat e_nuc) #read nao <<< $(cat num_ao) -#read nkpts <<< $(cat num_kpts) +#read nkpts <<< $(cat kpt_num) #read ndf <<< $(cat num_df) ##./create_ezfio_complex_4idx.py $ezfio $nel $natom $nmo $e_nucl $nao $nkpts ./create_ezfio_complex_3idx.py $ezfio $h5file #$nel $natom $nmo $e_nucl $nao $nkpts $ndf #Handle the orbital consitensy check qp_edit -c $ezfio &> /dev/null -cp $ezfio/{ao,mo}_basis/ao_md5 +#cp $ezfio/{ao,mo}_basis/ao_md5 #Read the integral -echo 'Read Integral' +#echo 'Read Integral' ################################################ diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index d76085d1..d909c99e 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -251,7 +251,7 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, with open('num_ao','w') as f: f.write(str(nao*Nk)) - with open('num_kpts','w') as f: + with open('kpt_num','w') as f: f.write(str(Nk)) # _ # |\ | _ | _ _. ._ |_) _ ._ | _ o _ ._ diff --git a/src/utils_complex/create_ezfio_complex_3idx.py b/src/utils_complex/create_ezfio_complex_3idx.py index 449050e8..a786f993 100755 --- a/src/utils_complex/create_ezfio_complex_3idx.py +++ b/src/utils_complex/create_ezfio_complex_3idx.py @@ -3,6 +3,7 @@ from ezfio import ezfio import h5py import sys +import numpy as np filename = sys.argv[1] h5filename = sys.argv[2] #num_elec, nucl_num, mo_num = map(int,sys.argv[2:5]) @@ -94,10 +95,6 @@ ezfio.set_ao_basis_ao_power(d) ezfio.set_ao_basis_ao_coef(d) ezfio.set_ao_basis_ao_expo(d) -#Dummy one -ao_md5 = '3b8b464dfc95f282129bde3efef3c502' -ezfio.set_ao_basis_ao_md5(ao_md5) -ezfio.set_mo_basis_ao_md5(ao_md5) ezfio.set_mo_basis_mo_num(mo_num) @@ -105,34 +102,63 @@ ezfio.set_mo_basis_mo_num(mo_num) #ezfio.set_mo_basis_mo_coef([ [0]*mo_num] * ao_num) ##ezfio.set_mo_basis_mo_coef_real(c_mo) +mo_coef_re0 = qph5['mo_basis/mo_coef_real'][()].T +mo_coef_im0 = qph5['mo_basis/mo_coef_imag'][()].T +mo_coef_cmplx0 = np.stack((mo_coef_re0,mo_coef_im0),axis=-1).tolist() -ezfio.set_mo_basis_mo_coef_real(qph5['mo_basis/mo_coef_real'][()].tolist()) -ezfio.set_mo_basis_mo_coef_imag(qph5['mo_basis/mo_coef_imag'][()].tolist()) +#ezfio.set_mo_basis_mo_coef_real(qph5['mo_basis/mo_coef_real'][()].tolist()) +#ezfio.set_mo_basis_mo_coef_imag(qph5['mo_basis/mo_coef_imag'][()].tolist()) +ezfio.set_mo_basis_mo_coef_complex(mo_coef_cmplx0) #maybe fix qp so we don't need this? ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num]) ezfio.set_nuclei_is_complex(True) +# fortran-ordered re,im parts +kin_ao_re0=qph5['ao_one_e_ints/ao_integrals_kinetic_real'][()].T +kin_ao_im0=qph5['ao_one_e_ints/ao_integrals_kinetic_imag'][()].T +#test where to stack? (axis=0 or -1?) +kin_ao_cmplx0=np.stack((kin_ao_re0,kin_ao_im0),axis=-1).tolist() -kin_ao_re=qph5['ao_one_e_ints/ao_integrals_kinetic_real'][()].T.tolist() -kin_ao_im=qph5['ao_one_e_ints/ao_integrals_kinetic_imag'][()].T.tolist() -ovlp_ao_re=qph5['ao_one_e_ints/ao_integrals_overlap_real'][()].T.tolist() -ovlp_ao_im=qph5['ao_one_e_ints/ao_integrals_overlap_imag'][()].T.tolist() -ne_ao_re=qph5['ao_one_e_ints/ao_integrals_n_e_real'][()].T.tolist() -ne_ao_im=qph5['ao_one_e_ints/ao_integrals_n_e_imag'][()].T.tolist() +ovlp_ao_re0=qph5['ao_one_e_ints/ao_integrals_overlap_real'][()].T +ovlp_ao_im0=qph5['ao_one_e_ints/ao_integrals_overlap_imag'][()].T +#test where to stack? (axis=0 or -1?) +ovlp_ao_cmplx0=np.stack((ovlp_ao_re0,ovlp_ao_im0),axis=-1).tolist() -ezfio.set_ao_one_e_ints_ao_integrals_kinetic(kin_ao_re) -ezfio.set_ao_one_e_ints_ao_integrals_kinetic_imag(kin_ao_im) -ezfio.set_ao_one_e_ints_ao_integrals_overlap(ovlp_ao_re) -ezfio.set_ao_one_e_ints_ao_integrals_overlap_imag(ovlp_ao_im) -ezfio.set_ao_one_e_ints_ao_integrals_n_e(ne_ao_re) -ezfio.set_ao_one_e_ints_ao_integrals_n_e_imag(ne_ao_im) +ne_ao_re0=qph5['ao_one_e_ints/ao_integrals_n_e_real'][()].T +ne_ao_im0=qph5['ao_one_e_ints/ao_integrals_n_e_imag'][()].T +#test where to stack? (axis=0 or -1?) +ne_ao_cmplx0=np.stack((ne_ao_re0,ne_ao_im0),axis=-1).tolist() + +kin_ao_re=kin_ao_re0.tolist() +kin_ao_im=kin_ao_im0.tolist() +ovlp_ao_re=ovlp_ao_re0.tolist() +ovlp_ao_im=ovlp_ao_im0.tolist() +ne_ao_re=ne_ao_re0.tolist() +ne_ao_im=ne_ao_im0.tolist() + +#kin_ao_c = np.stack(kin_ao_re0,kin_ao_im0 + +#ezfio.set_ao_one_e_ints_ao_integrals_kinetic(kin_ao_re) +#ezfio.set_ao_one_e_ints_ao_integrals_kinetic_imag(kin_ao_im) +ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_cmplx0) + +#ezfio.set_ao_one_e_ints_ao_integrals_overlap(ovlp_ao_re) +#ezfio.set_ao_one_e_ints_ao_integrals_overlap_imag(ovlp_ao_im) +ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_cmplx0) + +#ezfio.set_ao_one_e_ints_ao_integrals_n_e(ne_ao_re) +#ezfio.set_ao_one_e_ints_ao_integrals_n_e_imag(ne_ao_im) +ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_cmplx0) + +dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0)) +dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0)) +#ezfio.set_ao_two_e_ints_df_ao_integrals_real(dfao_re.tolist()) +#ezfio.set_ao_two_e_ints_df_ao_integrals_imag(dfao_im.tolist()) +dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist() +ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0) -dfao_re=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0)).tolist() -dfao_im=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0)).tolist() -ezfio.set_ao_two_e_ints_df_ao_integrals_real(dfao_re) -ezfio.set_ao_two_e_ints_df_ao_integrals_imag(dfao_im) #TODO: add check and only do this if ints exist #dfmo_re=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0)).tolist()