9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-02-16 05:34:04 +01:00

resolved dev kpts merge

This commit is contained in:
Kevin Gasperich 2020-06-16 10:38:27 -05:00
commit 464a6d70c4
152 changed files with 22640 additions and 1069 deletions

10
REPLACE
View File

@ -126,7 +126,6 @@ qp_name H_S2_u_0_bielec_nstates_openmp_work_3 -r H_S2_u_0_two_e_nstates_openmp_w
qp_name H_S2_u_0_bielec_nstates_openmp_work_4 -r H_S2_u_0_two_e_nstates_openmp_work_4
qp_name H_S2_u_0_bielec_nstates_openmp_work_$N_int
qp_name H_S2_u_0_bielec_nstates_openmp_work_$N_int -r "H_S2_u_0_two_e_nstates_openmp_work_$N_int"
qp_name H_S2_u_0_bielec_nstates_openmp_work_$N_int #-r "H_S2_u_0_two_e_nstates_openmp_work_$N_int"
qp_name H_S2_u_0_bielec_nstates_openmp_work -r H_S2_u_0_two_e_nstates_openmp_work
qp_name H_S2_u_0_bielec_nstates_openmp_work_ -r H_S2_u_0_two_e_nstates_openmp_work_
qp_name i_H_j_bielec -r i_H_j_two_e
@ -227,6 +226,7 @@ qp_name ref_bitmask_e_n_energy -r ref_bitmask_n_e_energy
qp_name read_ao_integrals_e_n -r read_ao_integrals_n_e
qp_name write_ao_integrals_e_n -r write_ao_integrals_n_e
qp_name psi_energy_bielec -r psi_energy_two_e
qp_name read_ao_integrals_e_n -r read_ao_integrals_n_e
qp_name read_ao_integrals --rename="read_ao_two_e_integrals"
qp_name read_ao_integrals --rename=read_ao_two_e_integrals
qp_name read_mo_integrals_erf -r read_mo_two_e_integrals_erf
@ -253,3 +253,11 @@ qp_name ezfio_set_mo_one_e_ints_mo_integrals_e_n -r ezfio_set_mo_one_e_ints_mo_i
qp_name ezfio_has_mo_one_e_ints_mo_integrals_e_n -r ezfio_has_mo_one_e_ints_mo_integrals_n_e
qp_name ezfio_has_mo_one_e_ints_io_mo_integrals_e_n -r ezfio_has_mo_one_e_ints_io_mo_integrals_n_e
qp_name ezfio_get_mo_one_e_ints_io_mo_integrals_e_n -r ezfio_get_mo_one_e_ints_io_mo_integrals_n_e
qp_name ao_ortho_canonical_coef_inv_complex -r ao_ortho_cano_coef_inv_cplx
qp_name fock_operator_closed_shell_ref_bitmask -r fock_op_cshell_ref_bitmask
qp_name fock_operator_closed_shell_ref_bitmask_complex -r fock_op_cshell_ref_bitmask_cplx
qp_name ao_ortho_canonical_coef_inv -r ao_ortho_cano_coef_inv
qp_name ao_ortho_cano_to_ao_complex -r ao_ortho_cano_to_ao_cplx
qp_name ao_ortho_lowdin_nucl_elec_integrals_complex -r ao_ortho_lowdin_n_e_ints_cplx
qp_name ao_ortho_canonical_nucl_elec_integrals_complex -r ao_ortho_cano_n_e_ints_cplx
qp_name ao_ortho_canonical_nucl_elec_integrals -r ao_ortho_cano_n_e_ints

1
configure vendored
View File

@ -365,7 +365,6 @@ EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file bse.tar.gz
pip install -e basis_set_exchange-*
EOF
elif [[ ${PACKAGE} = zlib ]] ; then
download ${ZLIB_URL} "${QP_ROOT}"/external/zlib.tar.gz

View File

@ -37,7 +37,9 @@ end = struct
} [@@deriving sexp]
;;
let get_default = Qpackage.get_ezfio_default "determinants";;
let get_default = Qpackage.get_ezfio_default "determinants"
let is_complex = lazy (Ezfio.get_nuclei_is_complex () )
let read_n_int () =
if not (Ezfio.has_determinants_n_int()) then
@ -48,12 +50,12 @@ end = struct
;
Ezfio.get_determinants_n_int ()
|> N_int_number.of_int
;;
let write_n_int n =
N_int_number.to_int n
|> Ezfio.set_determinants_n_int
;;
let read_bit_kind () =
@ -64,12 +66,12 @@ end = struct
;
Ezfio.get_determinants_bit_kind ()
|> Bit_kind.of_int
;;
let write_bit_kind b =
Bit_kind.to_int b
|> Ezfio.set_determinants_bit_kind
;;
let read_n_det () =
if not (Ezfio.has_determinants_n_det ()) then
@ -77,7 +79,7 @@ end = struct
;
Ezfio.get_determinants_n_det ()
|> Det_number.of_int
;;
let read_n_det_qp_edit () =
if not (Ezfio.has_determinants_n_det_qp_edit ()) then
@ -87,18 +89,18 @@ end = struct
end;
Ezfio.get_determinants_n_det_qp_edit ()
|> Det_number.of_int
;;
let write_n_det n =
Det_number.to_int n
|> Ezfio.set_determinants_n_det
;;
let write_n_det_qp_edit n =
let n_det = read_n_det () |> Det_number.to_int in
min n_det (Det_number.to_int n)
|> Ezfio.set_determinants_n_det_qp_edit
;;
let read_n_states () =
if not (Ezfio.has_determinants_n_states ()) then
@ -106,7 +108,7 @@ end = struct
;
Ezfio.get_determinants_n_states ()
|> States_number.of_int
;;
let write_n_states n =
let n_states =
@ -130,7 +132,7 @@ end = struct
Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_states |] ~data
|> Ezfio.set_determinants_state_average_weight
end
;;
let write_state_average_weight data =
let n_states =
@ -143,7 +145,7 @@ end = struct
in
Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_states |] ~data
|> Ezfio.set_determinants_state_average_weight
;;
let read_state_average_weight () =
let n_states =
@ -171,7 +173,7 @@ end = struct
|> Array.map Positive_float.of_float
in
(write_state_average_weight data; data)
;;
let read_expected_s2 () =
if not (Ezfio.has_determinants_expected_s2 ()) then
@ -186,12 +188,12 @@ end = struct
;
Ezfio.get_determinants_expected_s2 ()
|> Positive_float.of_float
;;
let write_expected_s2 s2 =
Positive_float.to_float s2
|> Ezfio.set_determinants_expected_s2
;;
let read_psi_coef ~read_only () =
if not (Ezfio.has_determinants_psi_coef ()) then
@ -200,19 +202,36 @@ end = struct
read_n_states ()
|> States_number.to_int
in
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| 1 ; n_states |]
~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. ))
(
if Lazy.force is_complex then
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| 1 ; n_states |]
~data:(List.init (2*n_states) (fun i -> if (i=0) then 1. else 0. ))
|> Ezfio.set_determinants_psi_coef
else
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| 2 ; 1 ; n_states |]
~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. ))
|> Ezfio.set_determinants_psi_coef_complex
)
end;
begin
if read_only then
Ezfio.get_determinants_psi_coef_qp_edit ()
begin
if Lazy.force is_complex then
Ezfio.get_determinants_psi_coef_complex_qp_edit ()
else
Ezfio.get_determinants_psi_coef_qp_edit ()
end
else
Ezfio.get_determinants_psi_coef ()
begin
if Lazy.force is_complex then
Ezfio.get_determinants_psi_coef_complex ()
else
Ezfio.get_determinants_psi_coef ()
end
end
|> Ezfio.flattened_ezfio
|> Array.map Det_coef.of_float
;;
let write_psi_coef ~n_det ~n_states c =
let n_det = Det_number.to_int n_det
@ -222,12 +241,23 @@ end = struct
and n_states =
States_number.to_int n_states
in
let r =
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
in
Ezfio.set_determinants_psi_coef r;
Ezfio.set_determinants_psi_coef_qp_edit r
;;
if Lazy.force is_complex then
begin
let r =
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| 2 ; n_det ; n_states |] ~data:c
in
Ezfio.set_determinants_psi_coef_complex r;
Ezfio.set_determinants_psi_coef_complex_qp_edit r
end
else
begin
let r =
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
in
Ezfio.set_determinants_psi_coef r;
Ezfio.set_determinants_psi_coef_qp_edit r
end
let read_psi_det ~read_only () =
@ -276,7 +306,7 @@ end = struct
|> Array.map (Determinant.of_int64_array
~n_int:(N_int_number.of_int n_int)
~alpha:n_alpha ~beta:n_beta )
;;
let write_psi_det ~n_int ~n_det d =
let data = Array.to_list d
@ -288,7 +318,7 @@ end = struct
in
Ezfio.set_determinants_psi_det r;
Ezfio.set_determinants_psi_det_qp_edit r
;;
let read ?(full=true) () =
@ -316,7 +346,7 @@ end = struct
else
(* No molecular orbitals, so no determinants *)
None
;;
let write ?(force=false)
{ n_int ;
@ -341,7 +371,7 @@ end = struct
write_psi_det ~n_int:n_int ~n_det:n_det psi_det
end;
write_state_average_weight state_average_weight
;;
let to_rst b =
@ -557,10 +587,8 @@ psi_det = %s
in
Generic_input_of_rst.evaluate_sexp t_of_sexp s
;;
let update_ndet n_det_new =
Printf.printf "Reducing n_det to %d\n" (Det_number.to_int n_det_new);
@ -596,7 +624,7 @@ psi_det = %s
{ det with n_det = (Det_number.of_int n_det_new) }
in
write ~force:true new_det
;;
let extract_state istate =
Printf.printf "Extracting state %d\n" (States_number.to_int istate);
@ -628,7 +656,7 @@ psi_det = %s
{ det with n_states = (States_number.of_int 1) }
in
write ~force:true new_det
;;
let extract_states range =
Printf.printf "Extracting states %s\n" (Range.to_string range);
@ -673,7 +701,7 @@ psi_det = %s
{ det with n_states = (States_number.of_int @@ List.length sorted_list) }
in
write ~force:true new_det
;;
end

View File

@ -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
@ -77,7 +73,7 @@ end = struct
let elec_alpha_num =
Ezfio.get_electrons_elec_alpha_num ()
in
let result =
let result =
Ezfio.get_mo_basis_mo_num ()
in
if result < elec_alpha_num then
@ -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.sub a (j*ao_num) (ao_num)
) )
else None
Array.init mo_num (fun j ->
Array.sub a (j*ao_num) (ao_num)
)
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
@ -262,7 +247,7 @@ mo_coef = %s
(b.mo_occ |> Array.to_list |> list_map
(MO_occ.to_string) |> String.concat ", " )
(b.mo_coef |> Array.map
(fun x-> Array.map MO_coef.to_string x |>
(fun x-> Array.map MO_coef.to_string x |>
Array.to_list |> String.concat "," ) |>
Array.to_list |> String.concat "\n" )
@ -300,40 +285,30 @@ 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
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_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
let write
let write
{ mo_num : MO_number.t ;
mo_label : MO_label.t;
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

View File

@ -166,6 +166,7 @@ let input_ezfio = "
let untouched = "
module MO_guess : sig
type t [@@deriving sexp]
val to_string : t -> string

View File

@ -67,3 +67,9 @@ doc: Use normalized primitive functions
interface: ezfio, provider
default: true
[ao_num_per_kpt]
type: integer
doc: Number of |AOs| per kpt
default: =(ao_basis.ao_num/nuclei.kpt_num)
interface: ezfio

View File

@ -0,0 +1,7 @@
BEGIN_PROVIDER [ integer, ao_num_per_kpt ]
implicit none
BEGIN_DOC
! number of aos per kpt.
END_DOC
ao_num_per_kpt = ao_num/kpt_num
END_PROVIDER

View File

@ -4,10 +4,16 @@ 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
[ao_integrals_n_e_kpts]
type: double precision
doc: Complex nucleus-electron integrals in |AO| basis set
size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_ao_integrals_n_e]
@ -23,10 +29,16 @@ 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
[ao_integrals_kinetic_kpts]
type: double precision
doc: Complex kinetic energy integrals in |AO| basis set
size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_ao_integrals_kinetic]
@ -42,10 +54,16 @@ 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
[ao_integrals_pseudo_kpts]
type: double precision
doc: Complex pseudopotential integrals in |AO| basis set
size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_ao_integrals_pseudo]
@ -61,10 +79,16 @@ 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
[ao_integrals_overlap_kpts]
type: double precision
doc: Complex overlap integrals in |AO| basis set
size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_ao_integrals_overlap]
@ -80,10 +104,16 @@ 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
[ao_one_e_integrals_kpts]
type: double precision
doc: Complex combined integrals in |AO| basis set
size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_ao_one_e_integrals]

View File

@ -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
@ -24,24 +27,85 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)]
!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
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
IF (read_ao_one_e_integrals) THEN
call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag)
call ezfio_get_ao_one_e_ints_ao_one_e_integrals_complex(ao_one_e_integrals_complex)
ELSE
print *, irp_here, ': Not yet implemented'
stop -1
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(ao_one_e_integrals_imag)
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
BEGIN_PROVIDER [ complex*16, ao_one_e_integrals_kpts,(ao_num_per_kpt,ao_num_per_kpt,kpt_num)]
&BEGIN_PROVIDER [ double precision, ao_one_e_integrals_diag_kpts,(ao_num_per_kpt,kpt_num)]
implicit none
integer :: j,k
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_kpts(ao_one_e_integrals_kpts)
else
ao_one_e_integrals_kpts = ao_integrals_n_e_kpts + ao_kinetic_integrals_kpts
if (do_pseudo) then
ao_one_e_integrals_kpts += ao_pseudo_integrals_kpts
endif
endif
do k = 1, kpt_num
do j = 1, ao_num_per_kpt
ao_one_e_integrals_diag_kpts(j,k) = dble(ao_one_e_integrals_kpts(j,j,k))
enddo
enddo
if (write_ao_one_e_integrals) then
call ezfio_set_ao_one_e_ints_ao_one_e_integrals_kpts(ao_one_e_integrals_kpts)
print *, 'AO one-e integrals written to disk'
endif
END_PROVIDER

View File

@ -84,13 +84,13 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_ortho_cano_coef_inv, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! ao_ortho_canonical_coef^(-1)
END_DOC
call get_inverse(ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),&
ao_num, ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1))
ao_num, ao_ortho_cano_coef_inv, size(ao_ortho_cano_coef_inv,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num,ao_num)]

View File

@ -0,0 +1,121 @@
!todo: add kpts
BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_complex, (ao_num,ao_cart_to_sphe_num) ]
implicit none
BEGIN_DOC
! complex version of ao_cart_to_sphe_coef
END_DOC
call zlacp2('A',ao_num,ao_cart_to_sphe_num, &
ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
ao_cart_to_sphe_coef_complex,size(ao_cart_to_sphe_coef_complex,1))
END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_overlap_complex, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ]
implicit none
BEGIN_DOC
! AO overlap matrix in the spherical basis set
END_DOC
complex*16, allocatable :: S(:,:)
allocate (S(ao_cart_to_sphe_num,ao_num))
call zgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, (1.d0,0.d0), &
ao_cart_to_sphe_coef_complex,size(ao_cart_to_sphe_coef_complex,1), &
ao_overlap_complex,size(ao_overlap_complex,1), (0.d0,0.d0), &
S, size(S,1))
call zgemm('N','N',ao_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, (1.d0,0.d0), &
S, size(S,1), &
ao_cart_to_sphe_coef_complex,size(ao_cart_to_sphe_coef_complex,1), (0.d0,0.d0), &
ao_cart_to_sphe_overlap_complex,size(ao_cart_to_sphe_overlap_complex,1))
deallocate(S)
END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_ortho_cano_coef_inv_cplx, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! ao_ortho_canonical_coef_complex^(-1)
END_DOC
call get_inverse_complex(ao_ortho_canonical_coef_complex,size(ao_ortho_canonical_coef_complex,1),&
ao_num, ao_ortho_cano_coef_inv_cplx, size(ao_ortho_cano_coef_inv_cplx,1))
END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_ortho_canonical_coef_complex, (ao_num,ao_num)]
&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num_complex ]
implicit none
BEGIN_DOC
! TODO: ao_ortho_canonical_num_complex should be the same as the real version
! maybe if the providers weren't linked we could avoid making a complex one?
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
END_DOC
integer :: i
ao_ortho_canonical_coef_complex = (0.d0,0.d0)
do i=1,ao_num
ao_ortho_canonical_coef_complex(i,i) = (1.d0,0.d0)
enddo
!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num)
!ao_ortho_canonical_num=ao_num
!return
if (ao_cartesian) then
ao_ortho_canonical_num_complex = ao_num
call ortho_canonical_complex(ao_overlap,size(ao_overlap,1), &
ao_num,ao_ortho_canonical_coef_complex,size(ao_ortho_canonical_coef_complex,1), &
ao_ortho_canonical_num_complex)
else
complex*16, allocatable :: S(:,:)
allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num))
S = (0.d0,0.d0)
do i=1,ao_cart_to_sphe_num
S(i,i) = (1.d0,0.d0)
enddo
ao_ortho_canonical_num_complex = ao_cart_to_sphe_num
call ortho_canonical_complex(ao_cart_to_sphe_overlap_complex, size(ao_cart_to_sphe_overlap_complex,1), &
ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num_complex)
call zgemm('N','N', ao_num, ao_ortho_canonical_num_complex, ao_cart_to_sphe_num, (1.d0,0.d0), &
ao_cart_to_sphe_coef_complex, size(ao_cart_to_sphe_coef_complex,1), &
S, size(S,1), &
(0.d0,0.d0), ao_ortho_canonical_coef_complex, size(ao_ortho_canonical_coef_complex,1))
deallocate(S)
endif
END_PROVIDER
BEGIN_PROVIDER [complex*16, ao_ortho_canonical_overlap_complex, (ao_ortho_canonical_num_complex,ao_ortho_canonical_num_complex)]
implicit none
BEGIN_DOC
! overlap matrix of the ao_ortho_canonical.
! Expected to be the Identity
END_DOC
integer :: i,j,k,l
complex*16 :: c
do j=1, ao_ortho_canonical_num_complex
do i=1, ao_ortho_canonical_num_complex
ao_ortho_canonical_overlap_complex(i,j) = (0.d0,0.d0)
enddo
enddo
do j=1, ao_ortho_canonical_num_complex
do k=1, ao_num
c = (0.d0,0.d0)
do l=1, ao_num
c += conjg(ao_ortho_canonical_coef_complex(l,j)) * ao_overlap_complex(l,k)
enddo
do i=1, ao_ortho_canonical_num_complex
ao_ortho_canonical_overlap_complex(i,j) += ao_ortho_canonical_coef_complex(k,i) * c
enddo
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,196 @@
!todo: add kpts
BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt)]
&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num_per_kpt ]
implicit none
BEGIN_DOC
! Coefficients to go from cartesian to spherical coordinates in the current
! basis set
END_DOC
integer :: i
integer, external :: ao_power_index
integer :: ibegin,j,k
integer :: prev
prev = 0
ao_cart_to_sphe_coef_kpts(:,:) = (0.d0,0.d0)
! Assume order provided by ao_power_index
i = 1
ao_cart_to_sphe_num_per_kpt = 0
do while (i <= ao_num_per_kpt)
select case ( ao_l(i) )
case (0)
ao_cart_to_sphe_num_per_kpt += 1
ao_cart_to_sphe_coef_kpts(i,ao_cart_to_sphe_num_per_kpt) = (1.d0,0.d0)
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef_kpts(i+j-1,ao_cart_to_sphe_num_per_kpt+k) = dcmplx(cart_to_sphe_$SHELL(j,k),0.d0)
enddo
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_num_per_kpt += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe_kpts : angular momentum too high'
end select
enddo
END_PROVIDER
!BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num_per_kpt ]
! implicit none
! ao_cart_to_sphe_num_per_kpt = ao_cart_to_sphe_num / kpt_num
!END_PROVIDER
!
!BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_kpts, (ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt) ]
! implicit none
! BEGIN_DOC
! ! complex version of ao_cart_to_sphe_coef for one k-point
! END_DOC
! call zlacp2('A',ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt, &
! ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
! ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1))
!END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_overlap_kpts, (ao_cart_to_sphe_num_per_kpt,ao_cart_to_sphe_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! AO overlap matrix in the spherical basis set
END_DOC
integer :: k
complex*16, allocatable :: S(:,:)
allocate (S(ao_cart_to_sphe_num_per_kpt,ao_num_per_kpt))
!todo: call with (:,:,k) vs (1,1,k)? is there a difference? does one create a temporary array?
do k=1, kpt_num
call zgemm('T','N',ao_cart_to_sphe_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt, (1.d0,0.d0), &
ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1), &
ao_overlap_kpts(:,:,k),size(ao_overlap_kpts,1), (0.d0,0.d0), &
S, size(S,1))
call zgemm('N','N',ao_cart_to_sphe_num_per_kpt,ao_cart_to_sphe_num_per_kpt,ao_num_per_kpt, (1.d0,0.d0), &
S, size(S,1), &
ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1), (0.d0,0.d0), &
ao_cart_to_sphe_overlap_kpts(:,:,k),size(ao_cart_to_sphe_overlap_kpts,1))
enddo
deallocate(S)
END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_ortho_cano_coef_inv_kpts, (ao_num_per_kpt,ao_num_per_kpt, kpt_num)]
implicit none
BEGIN_DOC
! ao_ortho_canonical_coef_complex^(-1)
END_DOC
integer :: k
do k=1, kpt_num
call get_inverse_complex(ao_ortho_canonical_coef_kpts,size(ao_ortho_canonical_coef_kpts,1),&
ao_num_per_kpt, ao_ortho_cano_coef_inv_kpts, size(ao_ortho_cano_coef_inv_kpts,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ complex*16, ao_ortho_canonical_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num)]
&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num_per_kpt, (kpt_num) ]
&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num_per_kpt_max ]
implicit none
BEGIN_DOC
! TODO: ao_ortho_canonical_num_complex should be the same as the real version
! maybe if the providers weren't linked we could avoid making a complex one?
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
END_DOC
integer :: i,k
ao_ortho_canonical_coef_kpts = (0.d0,0.d0)
do k=1,kpt_num
do i=1,ao_num
ao_ortho_canonical_coef_kpts(i,i,k) = (1.d0,0.d0)
enddo
enddo
!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num)
!ao_ortho_canonical_num=ao_num
!return
if (ao_cartesian) then
ao_ortho_canonical_num_per_kpt = ao_num_per_kpt
do k=1,kpt_num
call ortho_canonical_complex(ao_overlap_kpts(:,:,k),size(ao_overlap_kpts,1), &
ao_num_per_kpt,ao_ortho_canonical_coef_kpts(:,:,k),size(ao_ortho_canonical_coef_kpts,1), &
ao_ortho_canonical_num_per_kpt(k))
enddo
else
complex*16, allocatable :: S(:,:)
allocate(S(ao_cart_to_sphe_num_per_kpt,ao_cart_to_sphe_num_per_kpt))
do k=1,kpt_num
S = (0.d0,0.d0)
do i=1,ao_cart_to_sphe_num_per_kpt
S(i,i) = (1.d0,0.d0)
enddo
ao_ortho_canonical_num_per_kpt(k) = ao_cart_to_sphe_num_per_kpt
call ortho_canonical_complex(ao_cart_to_sphe_overlap_kpts, size(ao_cart_to_sphe_overlap_kpts,1), &
ao_cart_to_sphe_num_per_kpt, S, size(S,1), ao_ortho_canonical_num_per_kpt(k))
call zgemm('N','N', ao_num_per_kpt, ao_ortho_canonical_num_per_kpt(k), ao_cart_to_sphe_num_per_kpt, (1.d0,0.d0), &
ao_cart_to_sphe_coef_kpts, size(ao_cart_to_sphe_coef_kpts,1), &
S, size(S,1), &
(0.d0,0.d0), ao_ortho_canonical_coef_kpts(:,:,k), size(ao_ortho_canonical_coef_kpts,1))
enddo
deallocate(S)
endif
ao_ortho_canonical_num_per_kpt_max = maxval(ao_ortho_canonical_num_per_kpt)
END_PROVIDER
BEGIN_PROVIDER [complex*16, ao_ortho_canonical_overlap_kpts, (ao_ortho_canonical_num_per_kpt_max,ao_ortho_canonical_num_per_kpt_max,kpt_num)]
implicit none
BEGIN_DOC
! overlap matrix of the ao_ortho_canonical.
! Expected to be the Identity
END_DOC
integer :: i,j,k,l,kk
complex*16 :: c
do k=1,kpt_num
do j=1, ao_ortho_canonical_num_per_kpt_max
do i=1, ao_ortho_canonical_num_per_kpt_max
ao_ortho_canonical_overlap_kpts(i,j,k) = (0.d0,0.d0)
enddo
enddo
enddo
do kk=1,kpt_num
do j=1, ao_ortho_canonical_num_per_kpt(kk)
do k=1, ao_num_per_kpt
c = (0.d0,0.d0)
do l=1, ao_num_per_kpt
c += conjg(ao_ortho_canonical_coef_kpts(l,j,kk)) * ao_overlap_kpts(l,k,kk)
enddo
do i=1, ao_ortho_canonical_num_per_kpt(kk)
ao_ortho_canonical_overlap_kpts(i,j,kk) += ao_ortho_canonical_coef_kpts(k,i,kk) * c
enddo
enddo
enddo
enddo
END_PROVIDER

View File

@ -70,25 +70,55 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
implicit none
BEGIN_DOC
! Imaginary part of the overlap
END_DOC
ao_overlap_imag = 0.d0
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
BEGIN_PROVIDER [ complex*16, ao_overlap_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ]
implicit none
BEGIN_DOC
! Overlap for complex AOs
END_DOC
if (read_ao_integrals_overlap) then
call ezfio_get_ao_one_e_ints_ao_integrals_overlap_kpts(ao_overlap_kpts)
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_kpts(ao_overlap_kpts)
print *, 'AO overlap integrals written to disk'
endif
END_PROVIDER
@ -109,10 +139,15 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx
if (is_periodic) then
do j=1,ao_num
do i= 1,ao_num
ao_overlap_abs(i,j)= cdabs(ao_overlap_complex(i,j))
if (is_complex) then
ao_overlap_abs = 0.d0
integer :: k, ishift
do k=1,kpt_num
ishift = (k-1)*ao_num_per_kpt
do j=1,ao_num_per_kpt
do i= 1,ao_num_per_kpt
ao_overlap_abs(ishift+i,ishift+j)= cdabs(ao_overlap_kpts(i,j,k))
enddo
enddo
enddo
else
@ -175,6 +210,18 @@ BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]
ao_num,ao_num,S_inv_complex,size(S_inv_complex,1),lin_dep_cutoff)
END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_inv_kpts,(ao_num_per_kpt,ao_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! Inverse of the overlap matrix
END_DOC
integer :: k
do k=1,kpt_num
call get_pseudo_inverse_complex(ao_overlap_kpts(1,1,k), &
size(ao_overlap_kpts,1),ao_num_per_kpt,ao_num_per_kpt,S_inv_kpts(1,1,k),size(S_inv_kpts,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
BEGIN_DOC
@ -233,6 +280,125 @@ BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_half_inv_complex, (AO_num,AO_num) ]
BEGIN_DOC
! :math:`X = S^{-1/2}` obtained by SVD
END_DOC
implicit none
integer :: num_linear_dependencies
integer :: LDA, LDC
double precision, allocatable :: D(:)
complex*16, allocatable :: U(:,:),Vt(:,:)
integer :: info, i, j, k
double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6
LDA = size(AO_overlap_complex,1)
LDC = size(S_half_inv_complex,1)
allocate( &
U(LDC,AO_num), &
Vt(LDA,AO_num), &
D(AO_num))
call svd_complex( &
ao_overlap_complex,LDA, &
U,LDC, &
D, &
Vt,LDA, &
AO_num,AO_num)
num_linear_dependencies = 0
do i=1,AO_num
print*,D(i)
if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
D(i) = 0.d0
num_linear_dependencies += 1
else
ASSERT (D(i) > 0.d0)
D(i) = 1.d0/sqrt(D(i))
endif
do j=1,AO_num
S_half_inv_complex(j,i) = 0.d0
enddo
enddo
write(*,*) 'linear dependencies',num_linear_dependencies
do k=1,AO_num
if(D(k) /= 0.d0) then
do j=1,AO_num
do i=1,AO_num
S_half_inv_complex(i,j) = S_half_inv_complex(i,j) + U(i,k)*D(k)*Vt(k,j)
enddo
enddo
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_half_inv_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ]
BEGIN_DOC
! :math:`X = S^{-1/2}` obtained by SVD
END_DOC
implicit none
integer :: num_linear_dependencies
integer :: LDA, LDC
double precision, allocatable :: D(:)
complex*16, allocatable :: U(:,:),Vt(:,:)
integer :: info, i, j, k,kk
double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6
LDA = size(ao_overlap_kpts,1)
LDC = size(s_half_inv_kpts,1)
allocate( &
U(LDC,ao_num_per_kpt), &
Vt(LDA,ao_num_per_kpt), &
D(ao_num_per_kpt))
do kk=1,kpt_num
call svd_complex( &
ao_overlap_kpts(1,1,kk),LDA, &
U,LDC, &
D, &
Vt,LDA, &
ao_num_per_kpt,ao_num_per_kpt)
num_linear_dependencies = 0
do i=1,ao_num_per_kpt
print*,D(i)
if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
D(i) = 0.d0
num_linear_dependencies += 1
else
ASSERT (D(i) > 0.d0)
D(i) = 1.d0/sqrt(D(i))
endif
do j=1,ao_num_per_kpt
S_half_inv_kpts(j,i,kk) = 0.d0
enddo
enddo
write(*,*) 'linear dependencies, k: ',num_linear_dependencies,', ',kk
do k=1,ao_num_per_kpt
if(D(k) /= 0.d0) then
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
S_half_inv_kpts(i,j,kk) = S_half_inv_kpts(i,j,kk) + U(i,k)*D(k)*Vt(k,j)
enddo