diff --git a/REPLACE b/REPLACE index 04e50a6b..9485fb2a 100755 --- a/REPLACE +++ b/REPLACE @@ -1,5 +1,4 @@ # This file contains all the renamings that occured between qp1 and qp2. -# qp_name aa_operator_bielec -r aa_operator_two_e qp_name ac_operator_bielec -r ac_operator_two_e qp_name ao_bi_elec_integral_alpha -r ao_two_e_integral_alpha @@ -222,6 +221,11 @@ qp_name potential_sr_xc_beta_ao_lda --rename=potential_xc_beta_ao_sr_lda qp_name potential_sr_xc_beta_ao_lda --rename=potential_xc_beta_ao_sr_lda qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe +qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion +qp_name nucl_elec_ref_bitmask_energy -r ref_bitmask_n_e_energy +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 --rename="read_ao_two_e_integrals" qp_name read_ao_integrals --rename=read_ao_two_e_integrals @@ -240,3 +244,12 @@ qp_name write_ao_integrals --rename=write_ao_two_e_integrals qp_name write_mo_integrals_erf -r write_mo_two_e_integrals_erf qp_name write_mo_integrals --rename="write_mo_two_e_integrals" qp_name write_mo_integrals --rename=write_mo_two_e_integrals +qp_name io_mo_integrals_e_n -r io_mo_integrals_n_e +qp_name write_mo_integrals_e_n -r write_mo_integrals_n_e +qp_name read_mo_integrals_e_n -r read_mo_integrals_n_e +qp_name mo_integrals_e_n -r mo_integrals_n_e +qp_name ezfio_get_mo_one_e_ints_mo_integrals_e_n -r ezfio_get_mo_one_e_ints_mo_integrals_n_e +qp_name ezfio_set_mo_one_e_ints_mo_integrals_e_n -r ezfio_set_mo_one_e_ints_mo_integrals_n_e +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 diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index 12654aad..94435349 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -10,6 +10,7 @@ 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 @@ -24,6 +25,7 @@ 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" @@ -37,11 +39,18 @@ end = struct let reorder b ordering = - { b with mo_coef = - Array.map (fun mo -> - Array.init (Array.length mo) - (fun i -> mo.(ordering.(i))) - ) b.mo_coef + { b with + 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 ) } let read_ao_md5 () = @@ -121,6 +130,20 @@ end = struct 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 + let read () = if (Ezfio.has_mo_basis_mo_num ()) then @@ -130,6 +153,7 @@ 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 @@ -137,6 +161,7 @@ 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 = @@ -222,6 +247,7 @@ MO coefficients :: let to_string b = + (*TODO : add imaginary part here *) Printf.sprintf " mo_label = \"%s\" mo_num = %s @@ -244,12 +270,12 @@ mo_coef = %s let write_mo_num n = MO_number.to_int n |> Ezfio.set_mo_basis_mo_num - ;; + let write_mo_label a = MO_label.to_string a |> Ezfio.set_mo_basis_mo_label - ;; + let write_mo_class a = let mo_num = Array.length a in @@ -257,7 +283,7 @@ mo_coef = %s |> Array.to_list in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| mo_num |] ~data |> Ezfio.set_mo_basis_mo_class - ;; + let write_mo_occ a = let mo_num = Array.length a in @@ -265,12 +291,12 @@ mo_coef = %s |> Array.to_list in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| mo_num |] ~data |> Ezfio.set_mo_basis_mo_occ - ;; + let write_md5 a = MD5.to_string a |> Ezfio.set_mo_basis_ao_md5 - ;; + let write_mo_coef a = let mo_num = Array.length a in @@ -282,7 +308,24 @@ mo_coef = %s |> 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 + let write { mo_num : MO_number.t ; @@ -290,6 +333,7 @@ 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; @@ -297,8 +341,9 @@ 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 - ;; + end diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 930f7a77..9ef019fa 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -1,10 +1,16 @@ -[ao_integrals_e_n] +[ao_integrals_n_e] type: double precision doc: Nucleus-electron integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio -[io_ao_integrals_e_n] +[ao_integrals_n_e_imag] +type: double precision +doc: Imaginary part of the nucleus-electron integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_integrals_n_e] type: Disk_access doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml @@ -17,6 +23,12 @@ doc: Kinetic energy integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio +[ao_integrals_kinetic_imag] +type: double precision +doc: Imaginary part of the kinetic energy integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + [io_ao_integrals_kinetic] type: Disk_access doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ] @@ -30,6 +42,12 @@ doc: Pseudopotential integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio +[ao_integrals_pseudo_imag] +type: double precision +doc: Imaginary part of the pseudopotential integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + [io_ao_integrals_pseudo] type: Disk_access doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ] @@ -43,6 +61,12 @@ doc: Overlap integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio +[ao_integrals_overlap_imag] +type: double precision +doc: Imaginary part of the overlap integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + [io_ao_integrals_overlap] type: Disk_access doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ] @@ -56,6 +80,12 @@ doc: Combined integrals in |AO| basis set size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio +[ao_one_e_integrals_imag] +type: double precision +doc: Imaginary part of the combined integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + [io_ao_one_e_integrals] type: Disk_access doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ] 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 1df7fd37..c3084ae2 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 @@ -27,3 +27,24 @@ 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 + + IF (read_ao_one_e_integrals) THEN + call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag) + ELSE + print *, irp_here, ': Not yet implemented' + stop -1 + ENDIF + + IF (write_ao_one_e_integrals) THEN + call ezfio_set_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag) + 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 0339ca05..d7300936 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -70,6 +70,29 @@ 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 [ 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 +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] implicit none @@ -86,44 +109,52 @@ 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 - dim1=100 - lower_exp_val = 40.d0 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_x,overlap_y, overlap_z, overlap, & - !$OMP alpha, beta,i,j,dx) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - ao_overlap_abs(i,j)= 0.d0 - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) - call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1) - call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1) - ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z - enddo + 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)) + enddo enddo - enddo - enddo - !$OMP END PARALLEL DO + else + dim1=100 + lower_exp_val = 40.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B, & + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,& + !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + do i= 1,ao_num + ao_overlap_abs(i,j)= 0.d0 + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + do n = 1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(n,j) + do l = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(l,i) + call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) + call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1) + call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1) + ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif END_PROVIDER BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] @@ -134,6 +165,15 @@ BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) END_PROVIDER +BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Inverse of the overlap matrix + END_DOC + call get_pseudo_inverse_complex(ao_overlap_complex, & + size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1)) +END_PROVIDER + BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] BEGIN_DOC 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 b87e2a22..4f117deb 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -149,3 +149,25 @@ 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(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(ao_kinetic_integrals_imag) + 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 bb23f3b1..1d4cab7d 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -12,8 +12,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] 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_e_n) then - call ezfio_get_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e) + if (read_ao_integrals_n_e) then + call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e) print *, 'AO N-e integrals read from disk' else @@ -76,13 +76,36 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] !$OMP END DO !$OMP END PARALLEL endif - if (write_ao_integrals_e_n) then - call ezfio_set_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e) + if (write_ao_integrals_n_e) then + call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e) print *, 'AO N-e integrals written to disk' 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 [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)] BEGIN_DOC ! Nucleus-electron interaction in the |AO| basis set, per atom A. @@ -166,7 +189,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b double precision :: P_center(3) double precision :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi - double precision :: V_e_n,const_factor,dist_integral,tmp + double precision :: V_n_e,const_factor,dist_integral,tmp double precision :: accu,epsilo,rint integer :: n_pt_out,lmax include 'utils/constants.include.F' @@ -178,7 +201,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b (A_center(3)/=C_center(3))) then continue else - NAI_pol_mult = V_e_n(power_A(1),power_A(2),power_A(3), & + NAI_pol_mult = V_n_e(power_A(1),power_A(2),power_A(3), & power_B(1),power_B(2),power_B(3),alpha,beta) return endif @@ -476,7 +499,7 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim) endif end -double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta) +double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta) implicit none BEGIN_DOC ! Primitve nuclear attraction between the two primitves centered on the same atom. @@ -489,9 +512,9 @@ double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta) double precision :: alpha,beta double precision :: V_r, V_phi, V_theta if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then - V_e_n = 0.d0 + V_n_e = 0.d0 else - V_e_n = V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta) & + V_n_e = V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta) & * V_phi(a_x+b_x,a_y+b_y) & * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1) endif diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 9e729cd4..5272096d 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -21,7 +21,7 @@ subroutine two_e_integrals_index(i,j,k,l,i1) implicit none BEGIN_DOC ! Gives a unique index for i,j,k,l using permtuation symmetry. -! i <-> k, j <-> l, and (i,k) <-> (j,l) +! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems END_DOC integer, intent(in) :: i,j,k,l integer(key_kind), intent(out) :: i1 @@ -37,6 +37,8 @@ subroutine two_e_integrals_index(i,j,k,l,i1) i1 = i1+shiftr(i2*i2-i2,1) end + + subroutine two_e_integrals_index_reverse(i,j,k,l,i1) use map_module implicit none @@ -126,6 +128,155 @@ subroutine two_e_integrals_index_reverse(i,j,k,l,i1) end + + +subroutine ao_idx2_sq(i,j,ij) + implicit none + integer, intent(in) :: i,j + integer, intent(out) :: ij + if (ij) then + ij=(i-1)*(i-1)+2*j-mod(i,2) + else + ij=i*i + endif +end + +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 ao_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_2fold(i,j,k,l,i1) + use map_module + implicit none + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer :: ik,jl + + call ao_idx2_sq(i,k,ik) + call ao_idx2_sq(j,l,jl) + call ao_idx2_tri_key(ik,jl,i1) +end + +subroutine ao_idx2_sq_rev(i,k,ik) + BEGIN_DOC + ! reverse square compound index + END_DOC +! p = ceiling(dsqrt(dble(ik))) +! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1)))) +! if (mod(ik,2)==0) then +! k=p +! i=q +! else +! i=p +! k=q +! endif + integer, intent(in) :: ik + integer, intent(out) :: i,k + integer :: pq(0:1),i1,i2 + pq(0) = ceiling(dsqrt(dble(ik))) + pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1)))) + i1=mod(ik,2) + i2=mod(ik+1,2) + + k=pq(i1) + i=pq(i2) +end + +subroutine ao_idx2_tri_rev_key(i,k,ik) + use map_module + BEGIN_DOC + !return i<=k + END_DOC + integer(key_kind), intent(in) :: ik + integer, intent(out) :: i,k + integer(key_kind) :: tmp_k + k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) + tmp_k = k + i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1)) +end + +subroutine idx2_tri_rev_int(i,k,ik) + BEGIN_DOC + !return i<=k + END_DOC + integer, intent(in) :: ik + integer, intent(out) :: i,k + k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) + i = int(ik - ishft(k*k-k,-1)) +end + +subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1) + use map_module + implicit none + integer, intent(out) :: i(2),j(2),k(2),l(2) + integer(key_kind), intent(in) :: i1 + integer(key_kind) :: i0 + integer :: i2,i3 + i = 0 + call ao_idx2_tri_rev_key(i3,i2,i1) + + call ao_idx2_sq_rev(j(1),l(1),i2) + call ao_idx2_sq_rev(i(1),k(1),i3) + + !ijkl + i(2) = j(1) !jilk + j(2) = i(1) + k(2) = l(1) + l(2) = k(1) + +! i(3) = k(1) !klij complex conjugate +! j(3) = l(1) +! k(3) = i(1) +! l(3) = j(1) +! +! i(4) = l(1) !lkji complex conjugate +! j(4) = k(1) +! k(4) = j(1) +! l(4) = i(1) + + integer :: ii + if ( (i(1)==i(2)).and. & + (j(1)==j(2)).and. & + (k(1)==k(2)).and. & + (l(1)==l(2)) ) then + i(2) = 0 + endif +! This has been tested with up to 1000 AOs, and all the reverse indices are +! correct ! We can remove the test +! do ii=1,2 +! if (i(ii) /= 0) then +! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0) +! if (i1 /= i0) then +! print *, i1, i0 +! print *, i(ii), j(ii), k(ii), l(ii) +! stop 'two_e_integrals_index_reverse_2fold failed' +! endif +! endif +! enddo +end + + + + BEGIN_PROVIDER [ integer, ao_integrals_cache_min ] &BEGIN_PROVIDER [ integer, ao_integrals_cache_max ] implicit none @@ -144,28 +295,30 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ] END_DOC PROVIDE ao_two_e_integrals_in_map integer :: i,j,k,l,ii - integer(key_kind) :: idx + integer(key_kind) :: idx, idx2 real(integral_kind) :: integral - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) - do l=ao_integrals_cache_min,ao_integrals_cache_max - do k=ao_integrals_cache_min,ao_integrals_cache_max - do j=ao_integrals_cache_min,ao_integrals_cache_max - do i=ao_integrals_cache_min,ao_integrals_cache_max - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(ao_integrals_map,idx,integral) - ii = l-ao_integrals_cache_min - ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) - ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) - ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) - ao_integrals_cache(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO + real(integral_kind) :: tmp_re, tmp_im + integer(key_kind) :: idx_re,idx_im + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + do l=ao_integrals_cache_min,ao_integrals_cache_max + do k=ao_integrals_cache_min,ao_integrals_cache_max + do j=ao_integrals_cache_min,ao_integrals_cache_max + do i=ao_integrals_cache_min,ao_integrals_cache_max + !DIR$ FORCEINLINE + call two_e_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(ao_integrals_map,idx,integral) + ii = l-ao_integrals_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + ao_integrals_cache(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO END_PROVIDER @@ -207,6 +360,113 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result) result = tmp end +BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of AO integrals for fast access + END_DOC + PROVIDE ao_two_e_integrals_in_map + integer :: i,j,k,l,ii + integer(key_kind) :: idx1, idx2 + real(integral_kind) :: tmp_re, tmp_im + integer(key_kind) :: idx_re,idx_im + complex(integral_kind) :: integral + + + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral) + do l=ao_integrals_cache_min,ao_integrals_cache_max + do k=ao_integrals_cache_min,ao_integrals_cache_max + do j=ao_integrals_cache_min,ao_integrals_cache_max + do i=ao_integrals_cache_min,ao_integrals_cache_max + !DIR$ FORCEINLINE + call two_e_integrals_index_2fold(i,j,k,l,idx1) + !DIR$ FORCEINLINE + call two_e_integrals_index_2fold(k,l,i,j,idx2) + idx_re = min(idx1,idx2) + idx_im = max(idx1,idx2) + !DIR$ FORCEINLINE + call map_get(ao_integrals_map,idx_re,tmp_re) + if (idx_re /= idx_im) then + call map_get(ao_integrals_map,idx_im,tmp_im) + if (idx1 < idx2) then + integral = dcmplx(tmp_re,tmp_im) + else + integral = dcmplx(tmp_re,-tmp_im) + endif + else + tmp_im = 0.d0 + integral = dcmplx(tmp_re,tmp_im) + endif + + ii = l-ao_integrals_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + ao_integrals_cache_periodic(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result) + use map_module + implicit none + BEGIN_DOC + ! Gets one AO bi-electronic integral from the AO map + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx1,idx2 + real(integral_kind) :: tmp_re, tmp_im + integer(key_kind) :: idx_re,idx_im + type(map_type), intent(inout) :: map + integer :: ii + complex(integral_kind) :: tmp + PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min + !DIR$ FORCEINLINE + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then + tmp = (0.d0,0.d0) + else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then + tmp = (0.d0,0.d0) + else + ii = l-ao_integrals_cache_min + ii = ior(ii, k-ao_integrals_cache_min) + ii = ior(ii, j-ao_integrals_cache_min) + ii = ior(ii, i-ao_integrals_cache_min) + if (iand(ii, -64) /= 0) then + !DIR$ FORCEINLINE + call two_e_integrals_index_2fold(i,j,k,l,idx1) + !DIR$ FORCEINLINE + call two_e_integrals_index_2fold(k,l,i,j,idx2) + idx_re = min(idx1,idx2) + idx_im = max(idx1,idx2) + !DIR$ FORCEINLINE + call map_get(ao_integrals_map,idx_re,tmp_re) + if (idx_re /= idx_im) then + call map_get(ao_integrals_map,idx_im,tmp_im) + if (idx1 < idx2) then + tmp = dcmplx(tmp_re,tmp_im) + else + tmp = dcmplx(tmp_re,-tmp_im) + endif + else + tmp_im = 0.d0 + tmp = dcmplx(tmp_re,tmp_im) + endif + else + ii = l-ao_integrals_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + tmp = ao_integrals_cache_periodic(ii) + endif + result = tmp + endif +end + subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) use map_module @@ -237,6 +497,36 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) end + +subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) + use map_module + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All i are retrieved for j,k,l fixed. + ! physicist convention : + END_DOC + implicit none + integer, intent(in) :: j,k,l, sze + complex(integral_kind), intent(out) :: out_val(sze) + + integer :: i + integer(key_kind) :: hash + double precision :: thresh + PROVIDE ao_two_e_integrals_in_map ao_integrals_map + thresh = ao_integrals_threshold + + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + double precision :: get_ao_two_e_integral + do i=1,sze + out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) + enddo + +end + subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) use map_module implicit none @@ -407,81 +697,81 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) end -subroutine dump_ao_integrals(filename) - use map_module - implicit none - BEGIN_DOC - ! Save to disk the |AO| integrals - END_DOC - character*(*), intent(in) :: filename - integer(cache_key_kind), pointer :: key(:) - real(integral_kind), pointer :: val(:) - integer*8 :: i,j, n - if (.not.mpi_master) then - return - endif - call ezfio_set_work_empty(.False.) - open(unit=66,file=filename,FORM='unformatted') - write(66) integral_kind, key_kind - write(66) ao_integrals_map%sorted, ao_integrals_map%map_size, & - ao_integrals_map%n_elements - do i=0_8,ao_integrals_map%map_size - write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,& - ao_integrals_map%map(i)%n_elements - enddo - do i=0_8,ao_integrals_map%map_size - key => ao_integrals_map%map(i)%key - val => ao_integrals_map%map(i)%value - n = ao_integrals_map%map(i)%n_elements - write(66) (key(j), j=1,n), (val(j), j=1,n) - enddo - close(66) - -end +!subroutine dump_ao_integrals(filename) +! use map_module +! implicit none +! BEGIN_DOC +! ! Save to disk the |AO| integrals +! END_DOC +! character*(*), intent(in) :: filename +! integer(cache_key_kind), pointer :: key(:) +! real(integral_kind), pointer :: val(:) +! integer*8 :: i,j, n +! if (.not.mpi_master) then +! return +! endif +! call ezfio_set_work_empty(.False.) +! open(unit=66,file=filename,FORM='unformatted') +! write(66) integral_kind, key_kind +! write(66) ao_integrals_map%sorted, ao_integrals_map%map_size, & +! ao_integrals_map%n_elements +! do i=0_8,ao_integrals_map%map_size +! write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,& +! ao_integrals_map%map(i)%n_elements +! enddo +! do i=0_8,ao_integrals_map%map_size +! key => ao_integrals_map%map(i)%key +! val => ao_integrals_map%map(i)%value +! n = ao_integrals_map%map(i)%n_elements +! write(66) (key(j), j=1,n), (val(j), j=1,n) +! enddo +! close(66) +! +!end -integer function load_ao_integrals(filename) - implicit none - BEGIN_DOC - ! Read from disk the |AO| integrals - END_DOC - character*(*), intent(in) :: filename - integer*8 :: i - integer(cache_key_kind), pointer :: key(:) - real(integral_kind), pointer :: val(:) - integer :: iknd, kknd - integer*8 :: n, j - load_ao_integrals = 1 - open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') - read(66,err=98,end=98) iknd, kknd - if (iknd /= integral_kind) then - print *, 'Wrong integrals kind in file :', iknd - stop 1 - endif - if (kknd /= key_kind) then - print *, 'Wrong key kind in file :', kknd - stop 1 - endif - read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,& - ao_integrals_map%n_elements - do i=0_8, ao_integrals_map%map_size - read(66,err=99,end=99) ao_integrals_map%map(i)%sorted, & - ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements - call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size) - enddo - do i=0_8, ao_integrals_map%map_size - key => ao_integrals_map%map(i)%key - val => ao_integrals_map%map(i)%value - n = ao_integrals_map%map(i)%n_elements - read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) - enddo - call map_sort(ao_integrals_map) - load_ao_integrals = 0 - return - 99 continue - call map_deinit(ao_integrals_map) - 98 continue - stop 'Problem reading ao_integrals_map file in work/' - -end - +!integer function load_ao_integrals(filename) +! implicit none +! BEGIN_DOC +! ! Read from disk the |AO| integrals +! END_DOC +! character*(*), intent(in) :: filename +! integer*8 :: i +! integer(cache_key_kind), pointer :: key(:) +! real(integral_kind), pointer :: val(:) +! integer :: iknd, kknd +! integer*8 :: n, j +! load_ao_integrals = 1 +! open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') +! read(66,err=98,end=98) iknd, kknd +! if (iknd /= integral_kind) then +! print *, 'Wrong integrals kind in file :', iknd +! stop 1 +! endif +! if (kknd /= key_kind) then +! print *, 'Wrong key kind in file :', kknd +! stop 1 +! endif +! read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,& +! ao_integrals_map%n_elements +! do i=0_8, ao_integrals_map%map_size +! read(66,err=99,end=99) ao_integrals_map%map(i)%sorted, & +! ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements +! call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size) +! enddo +! do i=0_8, ao_integrals_map%map_size +! key => ao_integrals_map%map(i)%key +! val => ao_integrals_map%map(i)%value +! n = ao_integrals_map%map(i)%n_elements +! read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) +! enddo +! call map_sort(ao_integrals_map) +! load_ao_integrals = 0 +! return +! 99 continue +! call map_deinit(ao_integrals_map) +! 98 continue +! stop 'Problem reading ao_integrals_map file in work/' +! +!end +! 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 3ecdbbb2..a2bde897 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -354,10 +354,10 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals if (read_ao_two_e_integrals) then print*,'Reading the AO integrals' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - print*, 'AO integrals provided' - ao_two_e_integrals_in_map = .True. - return + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + print*, 'AO integrals provided' + ao_two_e_integrals_in_map = .True. + return endif print*, 'Providing the AO integrals' diff --git a/src/mo_basis/track_orb.irp.f b/src/bitmask/track_orb.irp.f similarity index 97% rename from src/mo_basis/track_orb.irp.f rename to src/bitmask/track_orb.irp.f index c49c545b..1cdde9cb 100644 --- a/src/mo_basis/track_orb.irp.f +++ b/src/bitmask/track_orb.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, mo_coef_begin_iteration, (ao_num,mo_num) ] BEGIN_DOC ! Void provider to store the coefficients of the |MO| basis at the beginning of the SCF iteration ! - ! Usefull to track some orbitals + ! Useful to track some orbitals END_DOC END_PROVIDER diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index a055aad3..874af46a 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -9,6 +9,12 @@ 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_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) + [mo_label] type: character*(64) doc: Label characterizing the MOS (Local, Canonical, Natural, *etc*) diff --git a/src/mo_basis/NEED b/src/mo_basis/NEED index 1a6a943c..5a504c20 100644 --- a/src/mo_basis/NEED +++ b/src/mo_basis/NEED @@ -1,3 +1,3 @@ ao_basis -ao_one_e_ints electrons +ao_one_e_ints diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index aa04fb01..73d33901 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -93,6 +93,59 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ] endif END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_coef_imag, (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(:,:) + logical :: exists + PROVIDE ezfio_filename + + + if (mpi_master) then + ! Coefs + call ezfio_has_mo_basis_mo_coef_imag(exists) + 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_imag with MPI' + endif + IRP_ENDIF + + if (exists) then + if (mpi_master) then + call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag) + write(*,*) 'Read mo_coef_imag' + endif + IRP_IF MPI + call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_imag with MPI' + endif + IRP_ENDIF + else + ! Orthonormalized AO basis + do i=1,mo_num + do j=1,ao_num + mo_coef_imag(j,i) = 0.d0 + enddo + enddo + endif +END_PROVIDER + BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ] implicit none BEGIN_DOC diff --git a/src/mo_one_e_ints/EZFIO.cfg b/src/mo_one_e_ints/EZFIO.cfg index 79a3b351..1cb77b6e 100644 --- a/src/mo_one_e_ints/EZFIO.cfg +++ b/src/mo_one_e_ints/EZFIO.cfg @@ -1,10 +1,10 @@ -[mo_integrals_e_n] +[mo_integrals_n_e] type: double precision doc: Nucleus-electron integrals in |MO| basis set size: (mo_basis.mo_num,mo_basis.mo_num) interface: ezfio -[io_mo_integrals_e_n] +[io_mo_integrals_n_e] type: Disk_access doc: Read/Write |MO| electron-nucleus attraction integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml diff --git a/src/mo_one_e_ints/pot_mo_ints.irp.f b/src/mo_one_e_ints/pot_mo_ints.irp.f index 90f7b06c..5cf4febd 100644 --- a/src/mo_one_e_ints/pot_mo_ints.irp.f +++ b/src/mo_one_e_ints/pot_mo_ints.irp.f @@ -4,8 +4,8 @@ BEGIN_PROVIDER [double precision, mo_integrals_n_e, (mo_num,mo_num)] ! Nucleus-electron interaction on the |MO| basis END_DOC - if (read_mo_integrals_e_n) then - call ezfio_get_mo_one_e_ints_mo_integrals_e_n(mo_integrals_n_e) + if (read_mo_integrals_n_e) then + call ezfio_get_mo_one_e_ints_mo_integrals_n_e(mo_integrals_n_e) print *, 'MO N-e integrals read from disk' else call ao_to_mo( & @@ -15,8 +15,8 @@ BEGIN_PROVIDER [double precision, mo_integrals_n_e, (mo_num,mo_num)] size(mo_integrals_n_e,1) & ) endif - if (write_mo_integrals_e_n) then - call ezfio_set_mo_one_e_ints_mo_integrals_e_n(mo_integrals_n_e) + if (write_mo_integrals_n_e) then + call ezfio_set_mo_one_e_ints_mo_integrals_n_e(mo_integrals_n_e) print *, 'MO N-e integrals written to disk' endif diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index ebdcfead..34c27c46 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -31,3 +31,9 @@ default: None doc: Nuclear repulsion (Computed automaticaly or Read in the |EZFIO|) type:double precision interface: ezfio + +[is_periodic] +type: logical +doc: If true, the calculation uses periodic boundary conditions +interface: ezfio, provider, ocaml +default: false diff --git a/src/tools/save_natorb.irp.f b/src/tools/save_natorb.irp.f index 88b28f06..f6331d13 100644 --- a/src/tools/save_natorb.irp.f +++ b/src/tools/save_natorb.irp.f @@ -19,7 +19,7 @@ program save_natorb call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None') call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None') call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None') - call ezfio_set_mo_one_e_ints_io_mo_integrals_e_n('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None') call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None') end diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index f0593d35..93b367aa 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -43,6 +43,690 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) end +subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) + implicit none + BEGIN_DOC + ! Compute A = U.D.Vt + ! + ! LDx : leftmost dimension of x + ! + ! Dimension of A is m x n + ! A,U,Vt are complex*16 + ! D is double precision + END_DOC + + integer, intent(in) :: LDA, LDU, LDVt, m, n + complex*16, intent(in) :: A(LDA,n) + complex*16, intent(out) :: U(LDU,m) + complex*16, intent(out) :: Vt(LDVt,n) + double precision,intent(out) :: D(min(m,n)) + complex*16,allocatable :: work(:) + double precision,allocatable :: rwork(:) + integer :: info, lwork, i, j, k, lrwork + + complex*16,allocatable :: A_tmp(:,:) + allocate (A_tmp(LDA,n)) + A_tmp = A + lrwork = 5*min(m,n) + + ! Find optimal size for temp arrays + allocate(work(1),rwork(lrwork)) + lwork = -1 + call zgesvd('A','A', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, rwork, info) + lwork = int(work(1)) + deallocate(work) + + allocate(work(lwork)) + call zgesvd('A','A', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, rwork, info) + deallocate(work,rwork,A_tmp) + + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + +end + +subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) + implicit none + BEGIN_DOC + ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. + ! + ! overlap : overlap matrix + ! + ! LDA : leftmost dimension of overlap array + ! + ! N : Overlap matrix is NxN (array is (LDA,N) ) + ! + ! C : Coefficients of the vectors to orthogonalize. On exit, + ! orthogonal vectors + ! + ! LDC : leftmost dimension of C + ! + ! m : Coefficients matrix is MxN, ( array is (LDC,N) ) + ! + END_DOC + + integer, intent(in) :: lda, ldc, n + integer, intent(out) :: m + complex*16, intent(in) :: overlap(lda,n) + complex*16, intent(inout) :: C(ldc,n) + complex*16, allocatable :: U(:,:) + complex*16, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + complex*16, allocatable :: S(:,:) + !DIR$ ATTRIBUTES ALIGN : 64 :: U, Vt, D + integer :: info, i, j + + if (n < 2) then + return + endif + + allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) + + call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) + + D(:) = dsqrt(D(:)) + m=n + do i=1,n + if ( D(i) >= 1.d-6 ) then + D(i) = 1.d0/D(i) + else + m = i-1 + print *, 'Removed Linear dependencies below:', 1.d0/D(m) + exit + endif + enddo + do i=m+1,n + D(i) = 0.d0 + enddo + + do i=1,m + if ( D(i) >= 1.d5 ) then + print *, 'Warning: Basis set may have linear dependence problems' + endif + enddo + + do j=1,n + do i=1,n + S(i,j) = U(i,j)*D(j) + enddo + enddo + + do j=1,n + do i=1,n + U(i,j) = C(i,j) + enddo + enddo + + call zgemm('N','N',n,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) + deallocate (U, Vt, D, S) + +end + + +subroutine ortho_qr_complex(A,LDA,m,n) + implicit none + BEGIN_DOC + ! Orthogonalization using Q.R factorization + ! + ! A : matrix to orthogonalize + ! + ! LDA : leftmost dimension of A + ! + ! n : Number of rows of A + ! + ! m : Number of columns of A + ! + END_DOC + integer, intent(in) :: m,n, LDA + complex*16, intent(inout) :: A(LDA,n) + + integer :: lwork, info + integer, allocatable :: jpvt(:) + complex*16, allocatable :: tau(:), work(:) + + allocate (jpvt(n), tau(n), work(1)) + LWORK=-1 + call zgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO ) + LWORK=2*int(WORK(1)) + deallocate(WORK) + allocate(WORK(LWORK)) + call zgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO ) + call zungqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO) + deallocate(WORK,jpvt,tau) +end + +subroutine ortho_qr_unblocked_complex(A,LDA,m,n) + implicit none + BEGIN_DOC + ! Orthogonalization using Q.R factorization + ! + ! A : matrix to orthogonalize + ! + ! LDA : leftmost dimension of A + ! + ! n : Number of rows of A + ! + ! m : Number of columns of A + ! + END_DOC + integer, intent(in) :: m,n, LDA + double precision, intent(inout) :: A(LDA,n) + + integer :: info + integer, allocatable :: jpvt(:) + double precision, allocatable :: tau(:), work(:) + + print *, irp_here, ': TO DO' + stop -1 + +! allocate (jpvt(n), tau(n), work(n)) +! call dgeqr2( m, n, A, LDA, TAU, WORK, INFO ) +! call dorg2r(m, n, n, A, LDA, tau, WORK, INFO) +! deallocate(WORK,jpvt,tau) +end + +subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) + implicit none + BEGIN_DOC + ! Compute C_new=C_old.S^-1/2 orthogonalization. + ! + ! overlap : overlap matrix + ! + ! LDA : leftmost dimension of overlap array + ! + ! N : Overlap matrix is NxN (array is (LDA,N) ) + ! + ! C : Coefficients of the vectors to orthogonalize. On exit, + ! orthogonal vectors + ! + ! LDC : leftmost dimension of C + ! + ! M : Coefficients matrix is MxN, ( array is (LDC,N) ) + ! + END_DOC + + integer, intent(in) :: LDA, ldc, n, m + complex*16, intent(in) :: overlap(lda,n) + complex*16, intent(inout) :: C(ldc,n) + complex*16, allocatable :: U(:,:) + complex*16, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + complex*16, allocatable :: S(:,:) + integer :: info, i, j, k + + if (n < 2) then + return + endif + + allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n)) + + call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(S,U,D,Vt,n,C,m) & + !$OMP PRIVATE(i,j,k) + + !$OMP DO + do i=1,n + if ( D(i) < 1.d-6 ) then + D(i) = 0.d0 + else + D(i) = 1.d0/dsqrt(D(i)) + endif + do j=1,n + S(j,i) = (0.d0,0.d0) + enddo + enddo + !$OMP END DO + + do k=1,n + if (D(k) /= 0.d0) then + !$OMP DO + do j=1,n + do i=1,n + S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + !$OMP END DO NOWAIT + endif + enddo + + !$OMP BARRIER + !$OMP DO + do j=1,n + do i=1,m + U(i,j) = C(i,j) + enddo + enddo + !$OMP END DO + + !$OMP END PARALLEL + + call zgemm('N','N',m,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) + + deallocate(U,Vt,S,D) +end + +subroutine get_inverse_complex(A,LDA,m,C,LDC) + implicit none + BEGIN_DOC + ! Returns the inverse of the square matrix A + END_DOC + integer, intent(in) :: m, LDA, LDC + complex*16, intent(in) :: A(LDA,m) + complex*16, intent(out) :: C(LDC,m) + + integer :: info,lwork + integer, allocatable :: ipiv(:) + complex*16,allocatable :: work(:) + allocate (ipiv(m), work(m*m)) + lwork = size(work) + C(1:m,1:m) = A(1:m,1:m) + call zgetrf(m,m,C,size(C,1),ipiv,info) + if (info /= 0) then + print *, info + stop 'error in inverse (zgetrf)' + endif + call zgetri(m,C,size(C,1),ipiv,work,lwork,info) + if (info /= 0) then + print *, info + stop 'error in inverse (zgetri)' + endif + deallocate(ipiv,work) +end + + +subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC) + implicit none + BEGIN_DOC + ! Find C = A^-1 + END_DOC + integer, intent(in) :: m,n, LDA, LDC + complex*16, intent(in) :: A(LDA,n) + complex*16, intent(out) :: C(LDC,m) + + double precision, allocatable :: D(:), rwork(:) + complex*16, allocatable :: U(:,:), Vt(:,:), work(:), A_tmp(:,:) + integer :: info, lwork + integer :: i,j,k + allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n),rwork(5*n)) + do j=1,n + do i=1,m + A_tmp(i,j) = A(i,j) + enddo + enddo + lwork = -1 + call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info) + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + lwork = int(real(work(1))) + deallocate(work) + allocate(work(lwork)) + call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info) + if (info /= 0) then + print *, info, ':: SVD failed' + stop 1 + endif + + do i=1,n + if (D(i)/D(1) > 1.d-10) then + D(i) = 1.d0/D(i) + else + D(i) = 0.d0 + endif + enddo + + C = (0.d0,0.d0) + do i=1,m + do j=1,n + do k=1,n + C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) + enddo + enddo + enddo + + deallocate(U,D,Vt,work,A_tmp,rwork) + +end + +subroutine lapack_diagd_diag_in_place_complex(eigvalues,eigvectors,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H(complex) + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax +! double precision, intent(out) :: eigvectors(nmax,n) + complex*16, intent(inout) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) +! double precision, intent(in) :: H(nmax,n) + complex*16,allocatable :: work(:) + integer ,allocatable :: iwork(:) +! complex*16,allocatable :: A(:,:) + double precision, allocatable :: rwork(:) + integer :: lrwork, lwork, info, i,j,l,k, liwork + +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + lwork = 2*n*n + 2*n + lrwork = 2*n*n + 5*n+ 1 + liwork = 5*n + 3 + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + + lwork = -1 + liwork = -1 + lrwork = -1 + ! get optimal work size + call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( real(work(1))) + liwork = iwork(1) + lrwork = int(rwork(1)) + deallocate (work,iwork,rwork) + + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + deallocate(work,iwork,rwork) + + + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'ZHEEVD Failed; calling ZHEEV' + lwork = 2*n - 1 + lrwork = 3*n - 2 + allocate(work(lwork),rwork(lrwork)) + lwork = -1 + call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info) + if (info < 0) then + print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info) + if (info /= 0 ) then + write(*,*)'ZHEEV Failed' + stop 1 + endif + deallocate(work,rwork) + end if + +end + +subroutine lapack_diagd_diag_complex(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H(complex) + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax +! double precision, intent(out) :: eigvectors(nmax,n) + complex*16, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) +! double precision, intent(in) :: H(nmax,n) + complex*16, intent(in) :: H(nmax,n) + double precision, allocatable :: eigenvalues(:) + complex*16,allocatable :: work(:) + integer ,allocatable :: iwork(:) + complex*16,allocatable :: A(:,:) + double precision, allocatable :: rwork(:) + integer :: lrwork, lwork, info, i,j,l,k, liwork + + allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + A=H + lwork = 2*n*n + 2*n + lrwork = 2*n*n + 5*n+ 1 + liwork = 5*n + 3 + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + + lwork = -1 + liwork = -1 + lrwork = -1 + ! get optimal work size + call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( real(work(1))) + liwork = iwork(1) + lrwork = int(rwork(1)) + deallocate (work,iwork,rwork) + + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + deallocate(work,iwork,rwork) + + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'ZHEEVD Failed; calling ZHEEV' + lwork = 2*n - 1 + lrwork = 3*n - 2 + allocate(work(lwork),rwork(lrwork)) + lwork = -1 + call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info) + if (info < 0) then + print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info) + if (info /= 0 ) then + write(*,*)'ZHEEV Failed' + stop 1 + endif + deallocate(work,rwork) + end if + + eigvectors = (0.d0,0.d0) + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + +subroutine lapack_diagd_complex(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H(complex) + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax +! double precision, intent(out) :: eigvectors(nmax,n) + complex*16, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) +! double precision, intent(in) :: H(nmax,n) + complex*16, intent(in) :: H(nmax,n) + double precision, allocatable :: eigenvalues(:) + complex*16,allocatable :: work(:) + integer ,allocatable :: iwork(:) + complex*16,allocatable :: A(:,:) + double precision, allocatable :: rwork(:) + integer :: lrwork, lwork, info, i,j,l,k, liwork + + allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + A=H + lwork = 2*n*n + 2*n + lrwork = 2*n*n + 5*n+ 1 + liwork = 5*n + 3 + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + + lwork = -1 + liwork = -1 + lrwork = -1 + call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = max(int( work( 1 ) ),lwork) + liwork = iwork(1) + lrwork = max(int(rwork(1),4),lrwork) + deallocate (work,iwork,rwork) + + allocate (work(lwork),iwork(liwork),rwork(lrwork)) + call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, lrwork, iwork, liwork, info ) + deallocate(work,iwork,rwork) + + + if (info < 0) then + print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'ZHEEVD Failed' + stop 1 + end if + + eigvectors = (0.d0,0.d0) + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + +subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H (complex) + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax + complex*16, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) + complex*16, intent(in) :: H(nmax,n) + double precision,allocatable :: eigenvalues(:) + complex*16,allocatable :: work(:) + complex*16,allocatable :: A(:,:) + double precision,allocatable :: rwork(:) + integer :: lwork, info, i,j,l,k,lrwork + + allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + A=H + !lwork = 2*n*n + 6*n+ 1 + lwork = 2*n - 1 + lrwork = 3*n - 2 + allocate (work(lwork),rwork(lrwork)) + + lwork = -1 + call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, info ) + if (info < 0) then + print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( work( 1 ) ) + deallocate (work) + + allocate (work(lwork)) + call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + rwork, info ) + deallocate(work,rwork) + + if (info < 0) then + print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'ZHEEV Failed : ', info + do i=1,n + do j=1,n + print *, H(i,j) + enddo + enddo + stop 1 + end if + + eigvectors = (0.d0,0.d0) + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + +subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda) + implicit none + BEGIN_DOC +! performs u1 += u0 * matrix + END_DOC + integer, intent(in) :: sze,lda + complex*16, intent(in) :: u0(sze) + complex*16, intent(inout) :: u1(sze) + complex*16, intent(in) :: matrix(lda,sze) + integer :: i,j + integer :: incx,incy + incx = 1 + incy = 1 + !call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy) + call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy) +end + subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) implicit none BEGIN_DOC @@ -356,6 +1040,8 @@ subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC) end + + subroutine find_rotation(A,LDA,B,m,C,n) implicit none BEGIN_DOC diff --git a/src/utils_periodic/NEED b/src/utils_periodic/NEED new file mode 100644 index 00000000..49c1d899 --- /dev/null +++ b/src/utils_periodic/NEED @@ -0,0 +1,4 @@ +ao_two_e_ints +ao_one_e_ints +mo_two_e_ints + diff --git a/src/utils_periodic/import_integrals_ao_periodic.irp.f b/src/utils_periodic/import_integrals_ao_periodic.irp.f new file mode 100644 index 00000000..79eb8fe0 --- /dev/null +++ b/src/utils_periodic/import_integrals_ao_periodic.irp.f @@ -0,0 +1,137 @@ +program print_integrals + print *, 'Number of AOs?' + read(*,*) ao_num + TOUCH ao_num + call run +end + +subroutine run + use map_module + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j,k,l + double precision :: integral + double precision, allocatable :: A(:,:), B(:,:) + double precision :: tmp_re, tmp_im + + integer :: n_integrals + integer(key_kind), allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_values(:) + + call ezfio_set_ao_basis_ao_num(ao_num) + + allocate (A(ao_num,ao_num), B(ao_num,ao_num) ) + + A(1,1) = huge(1.d0) + iunit = getunitandopen('E.qp','r') + read (iunit,*,end=9) A(1,1) + 9 continue + close(iunit) + if (A(1,1) /= huge(1.d0)) then + call ezfio_set_nuclei_nuclear_repulsion(A(1,1)) + call ezfio_set_nuclei_io_nuclear_repulsion("Read") + endif + + A = 0.d0 + B = 0.d0 + iunit = getunitandopen('T.qp','r') + do + read (iunit,*,end=10) i,j, tmp_re, tmp_im + A(i,j) = tmp_re + B(i,j) = tmp_im + if (i.ne.j) then + A(j,i) = tmp_re + B(j,i) = -tmp_im + endif + enddo + 10 continue + close(iunit) + call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(B(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read") + + A = 0.d0 + B = 0.d0 + iunit = getunitandopen('S.qp','r') + do + read (iunit,*,end=11) i,j, tmp_re, tmp_im + A(i,j) = tmp_re + B(i,j) = tmp_im + if (i.ne.j) then + A(j,i) = tmp_re + B(j,i) = -tmp_im + endif + enddo + 11 continue + close(iunit) + call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(B(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read") + + A = 0.d0 + B = 0.d0 + iunit = getunitandopen('P.qp','r') + do + read (iunit,*,end=14) i,j, tmp_re, tmp_im + A(i,j) = tmp_re + B(i,j) = tmp_im + if (i.ne.j) then + A(j,i) = tmp_re + B(j,i) = -tmp_im + endif + enddo + 14 continue + close(iunit) + call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num)) + call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(B(1:ao_num,1:ao_num)) + call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read") + + A = 0.d0 + B = 0.d0 + iunit = getunitandopen('V.qp','r') + do + read (iunit,*,end=12) i,j, tmp_re, tmp_im + A(i,j) = tmp_re + B(i,j) = tmp_im + if (i.ne.j) then + A(j,i) = tmp_re + B(j,i) = -tmp_im + endif + enddo + 12 continue + close(iunit) + call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num)) + call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read") + +! allocate(buffer_i(ao_num**3), buffer_values(ao_num**3)) +! iunit = getunitandopen('W.qp','r') +! n_integrals=0 +! buffer_values = 0.d0 +! do +! read (iunit,*,end=13) i,j,k,l, integral +! n_integrals += 1 +! call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) ) +! buffer_values(n_integrals) = integral +! if (n_integrals == size(buffer_i)) then +! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values) +! n_integrals = 0 +! endif +! enddo +! 13 continue +! close(iunit) +! +! if (n_integrals > 0) then +! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values) +! endif +! +! call map_sort(ao_integrals_map) +! call map_unique(ao_integrals_map) +! +! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) +! call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + +end diff --git a/src/utils_periodic/import_mo_coef_periodic.irp.f b/src/utils_periodic/import_mo_coef_periodic.irp.f new file mode 100644 index 00000000..bd41f776 --- /dev/null +++ b/src/utils_periodic/import_mo_coef_periodic.irp.f @@ -0,0 +1,27 @@ +program import_mo_coef_periodic + + PROVIDE ezfio_filename + call run +end + +subroutine run + use map_module + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j + double precision :: int_re, int_im + + + iunit = getunitandopen('C.qp','r') + do + read (iunit,*,end=10) i,j, mo_coef(i,j), mo_coef_imag(i,j) + enddo + 10 continue + close(iunit) + mo_label = "None" + call save_mos + +end