From 9d3900c7ee39df5fbfc16ca955589ea02b479611 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 27 Jan 2016 17:15:57 +0100 Subject: [PATCH 1/3] qp_module is now case insensitive --- configure | 2 +- scripts/module/qp_module.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/configure b/configure index a6c5bc2e..2344264a 100755 --- a/configure +++ b/configure @@ -533,7 +533,7 @@ def recommendation(): print " source {0}".format(path) print "" print "Then, install the modules you want to install using :" - print " qp_install_module.py " + print " qp_module.py" print "" print "Finally :" print " ninja" diff --git a/scripts/module/qp_module.py b/scripts/module/qp_module.py index 3eda69f3..63649df5 100755 --- a/scripts/module/qp_module.py +++ b/scripts/module/qp_module.py @@ -175,7 +175,11 @@ def main(arguments): d_child = d_local.copy() d_child.update(d_plugin) - l_name = arguments[""] + normalize_case = {} + for name in d_local.keys() + d_plugin.keys(): + normalize_case [ name.lower() ] = name + + l_name = [ normalize_case[name.lower()] for name in arguments[""] ] for name in l_name: if name in d_local: @@ -206,6 +210,7 @@ def main(arguments): raise print "[ OK ]" + print "" print "You can now compile as usual" print "`cd {0} ; ninja` for exemple".format(QP_ROOT) print " or --in developement mode-- you can cd in a directory and compile here" From ccc061fbeb0099f4abacbf0ee3836af6727cb66f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 28 Jan 2016 14:50:24 +0100 Subject: [PATCH 2/3] Introduced MR-CEPA. Not working yet --- plugins/MRCC_Utils/H_apply.irp.f | 22 ++ plugins/MRCC_Utils/mrcc_utils.irp.f | 4 + plugins/MRCC_Utils/mrcepa_dress.irp.f | 260 ++++++++++++++++++++++++ plugins/MRCC_Utils/mrcepa_general.irp.f | 97 +++++++++ 4 files changed, 383 insertions(+) create mode 100644 plugins/MRCC_Utils/mrcepa_dress.irp.f create mode 100644 plugins/MRCC_Utils/mrcepa_general.irp.f diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 7314674c..1cafc8de 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -24,5 +24,27 @@ s.data["size_max"] = "3072" print s +s = H_apply("mrcepa") +s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" +s.data["declarations"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["decls_main"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["finalization"] = "" +s.data["copy_buffer"] = "" +s.data["generate_psi_guess"] = "" +s.data["size_max"] = "3072" +# print s + + END_SHELL diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 8fe6c411..1e2f974d 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -26,7 +26,11 @@ ! TODO --- Test perturbatif ------ do k=1,N_states lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + ! TODO : i_h_psi peut sortir de la boucle? call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif tmp = psi_non_ref_coef(i,k)/ihpsi_current(k) i_pert = 0 ! Perturbation only if 1st order < 0.5 x second order diff --git a/plugins/MRCC_Utils/mrcepa_dress.irp.f b/plugins/MRCC_Utils/mrcepa_dress.irp.f new file mode 100644 index 00000000..9789b818 --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_dress.irp.f @@ -0,0 +1,260 @@ +use omp_lib +use bitmasks + +subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_ref, Ndet_non_ref + double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(inout) :: delta_ii_(Ndet_ref,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l + integer :: degree_alpha(psi_det_size) + integer :: idx_alpha(0:psi_det_size) + logical :: good, fullMatch + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref ,degree + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hia(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer(bit_kind) :: tmp_det_0(Nint,2) + integer :: iint, ipos + integer :: i_state, i_sd, k_sd, l_sd, i_I, i_alpha + + integer(bit_kind),allocatable :: miniList(:,:,:) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer,allocatable :: idx_miniList(:) + integer :: N_miniList, ni, leng + integer(bit_kind) :: isum + + double precision :: hia + integer, allocatable :: index_sorted(:) + + + leng = max(N_det_generators, N_det_non_ref) + allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det)) + + !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) + + if(fullMatch) then + return + end if + + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + + allocate (dIa_hia(N_states,Ndet_non_ref)) + + ! |I> + + ! |alpha> + + if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint) + end if + + + do i_alpha=1,N_tq + ! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + + integer, external :: get_index_in_psi_det_sorted_bit + index_sorted = huge(-1) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j) + end do + + ! |I> + do i_I=1,N_det_ref + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint) + if (degree > 4) then + cycle + endif + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + + !TODO: MR + do i_sd=1,idx_alpha(0) + call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint) + if (degree > 2) then + cycle + endif + call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det_0 = 0_bit_kind + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + endif + + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia) + + ! |alpha> + do k_sd=1,idx_alpha(0) + call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) + if (degree > 2) then + cycle + endif + + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det = 0_bit_kind + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + isum = 0_bit_kind + do iint = 1,N_int + isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) & + + iand(tmp_det(iint,2), tmp_det_0(iint,2)) + enddo + + if (isum /= 0_bit_kind) then + cycle + endif + + ! + ! + call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) + do i_state=1,N_states + dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + enddo + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do k=1,N_int + tmp_det(k,1) = psi_ref(k,1,i_I) + tmp_det(k,2) = psi_ref(k,2,i_I) + enddo + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + + +! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) ) +! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl) +! do i_state=1,N_states +! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2 +! enddo + + do l_sd=1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) + if (degree == 0) then + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + enddo + exit + endif + enddo + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,N_states + ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state) + enddo + + k_sd = idx_alpha(i_sd) + do i_state=1,N_states + dIa_hia(i_state,k_sd) = dIa(i_state) * hia + enddo + + call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states + delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd) + + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) + else + delta_ii_(i_I,i_state) = 0.d0 + endif + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) + enddo + enddo + + enddo + deallocate (dIa_hia,index_sorted) + deallocate(miniList, idx_miniList) +end + + + + + + + + diff --git a/plugins/MRCC_Utils/mrcepa_general.irp.f b/plugins/MRCC_Utils/mrcepa_general.irp.f new file mode 100644 index 00000000..3479548b --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_general.irp.f @@ -0,0 +1,97 @@ +subroutine run_mrcepa + implicit none + call set_generators_bitmasks_as_holes_and_particles + call mrcepa_iterations +end + +subroutine mrcepa_iterations + implicit none + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration,i_oscillations + double precision :: E_past(4) + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + j = 1 + i_oscillations = 0 + do while (delta_E > 1.d-7) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + + E_past(j) = E_new + j +=1 + if(j>4)then + j=1 + endif + if(iteration > 4) then + if(delta_E > 1.d-10)then + if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then + print*,'OSCILLATIONS !!!' + oscillations = .True. + i_oscillations +=1 + lambda_mrcc_tmp = lambda_mrcc + endif + endif + endif + call save_wavefunction +! if (i_oscillations > 5) then +! exit +! endif + if (iteration > 200) then + exit + endif + print*,'------------' + print*,'VECTOR' + do i = 1, N_det_ref + print*,'' + print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) + print*,'delta_ii(i,1) = ',delta_ii(i,1) + enddo + print*,'------------' + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + +end + +subroutine set_generators_bitmasks_as_holes_and_particles + implicit none + integer :: i,k + do k = 1, N_generators_bitmask + do i = 1, N_int + ! Pure single part + generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha + generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha + generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta + generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta + + ! Double excitation + generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha + generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha + generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta + generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta + + generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha + generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha + generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta + generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta + + enddo + enddo + touch generators_bitmask + + + +end From d94138cfedeede81283819a3f65de71ed9bd41ba Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Feb 2016 00:37:03 +0100 Subject: [PATCH 3/3] Corrected bug when writing pseudo to EZFIO --- ocaml/Qputils.ml | 20 +- ocaml/qp_create_ezfio_from_xyz.ml | 553 ++++++++++++++++-------------- 2 files changed, 308 insertions(+), 265 deletions(-) diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 7f840fde..36dc1102 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core.Std (* let rec transpose = function @@ -24,5 +24,21 @@ let input_to_sexp s = print_endline ("("^result^")"); "("^result^")" |> Sexp.of_string -;; + +let rmdir dirname = + let rec remove_one dir = + Sys.chdir dir; + Sys.readdir "." + |> Array.iter ~f:(fun x -> + match (Sys.is_directory x, Sys.is_file x) with + | (`Yes, _) -> remove_one x + | (_, `Yes) -> Sys.remove x + | _ -> failwith ("Unable to remove file "^x^".") + ); + Sys.chdir ".."; + Unix.rmdir dir + in + remove_one dirname + + diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index bb7fd847..a6820d88 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -1,6 +1,6 @@ -open Qputils;; -open Qptypes;; -open Core.Std;; +open Qputils +open Qptypes +open Core.Std let spec = let open Command.Spec in @@ -316,289 +316,316 @@ let run ?o b c d m p cart xyz_file = if Sys.file_exists_exn ezfio_file then failwith (ezfio_file^" already exists"); - (* Create EZFIO *) - Ezfio.set_file ezfio_file; + let write_file () = + (* Create EZFIO *) + Ezfio.set_file ezfio_file; - (* Write Pseudo *) - let pseudo = - List.map nuclei ~f:(fun x -> - match pseudo_channel x.Atom.element with - | Some channel -> Pseudo.read_element channel x.Atom.element - | None -> Pseudo.empty x.Atom.element - ) - in + (* Write Pseudo *) + let pseudo = + List.map nuclei ~f:(fun x -> + match pseudo_channel x.Atom.element with + | Some channel -> Pseudo.read_element channel x.Atom.element + | None -> Pseudo.empty x.Atom.element + ) + in - let molecule = - let n_elec_to_remove = - List.fold pseudo ~init:0 ~f:(fun accu x -> - accu + (Positive_int.to_int x.Pseudo.n_elec)) - in - { Molecule.elec_alpha = - (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) - - n_elec_to_remove/2 - |> Elec_alpha_number.of_int; - Molecule.elec_beta = - (Elec_beta_number.to_int molecule.Molecule.elec_beta) - - (n_elec_to_remove - n_elec_to_remove/2) - |> Elec_beta_number.of_int; - Molecule.nuclei = - let charges = - List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec - |> Float.of_int) - |> Array.of_list + let molecule = + let n_elec_to_remove = + List.fold pseudo ~init:0 ~f:(fun accu x -> + accu + (Positive_int.to_int x.Pseudo.n_elec)) in - List.mapi molecule.Molecule.nuclei ~f:(fun i x -> - { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) - |> Charge.of_float } - ) - } - in - let nuclei = - molecule.Molecule.nuclei @ dummy - in - - - (* Write Electrons *) - Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int - molecule.Molecule.elec_alpha ) ; - Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int - molecule.Molecule.elec_beta ) ; - - (* Write Nuclei *) - let labels = - List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei - and charges = - List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei - and coords = - (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in - let nucl_num = (List.length labels) in - Ezfio.set_nuclei_nucl_num nucl_num ; - Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:labels); - Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:charges); - Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); - - - (* Write pseudopotential *) - let () = - match p with - | None -> Ezfio.set_pseudo_do_pseudo false - | _ -> Ezfio.set_pseudo_do_pseudo true - in - - let klocmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.local + { Molecule.elec_alpha = + (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) + - n_elec_to_remove/2 + |> Elec_alpha_number.of_int; + Molecule.elec_beta = + (Elec_beta_number.to_int molecule.Molecule.elec_beta) + - (n_elec_to_remove - n_elec_to_remove/2) + |> Elec_beta_number.of_int; + Molecule.nuclei = + let charges = + List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec + |> Float.of_int) + |> Array.of_list + in + List.mapi molecule.Molecule.nuclei ~f:(fun i x -> + { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) + |> Charge.of_float } + ) + } in - if (x > accu) then x - else accu - ) - and kmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.non_local + let nuclei = + molecule.Molecule.nuclei @ dummy in - if (x > accu) then x - else accu - ) - and lmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + + + (* Write Electrons *) + Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int + molecule.Molecule.elec_alpha ) ; + Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int + molecule.Molecule.elec_beta ) ; + + (* Write Nuclei *) + let labels = + List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei + and charges = + List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei + and coords = + (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in + let nucl_num = (List.length labels) in + Ezfio.set_nuclei_nucl_num nucl_num ; + Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:labels); + Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:charges); + Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); + + + (* Write pseudopotential *) + let () = + match p with + | None -> Ezfio.set_pseudo_do_pseudo false + | _ -> Ezfio.set_pseudo_do_pseudo true + in + + let klocmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> let x = - Positive_int.to_int x.Pseudo.Primitive_non_local.proj + List.length x.Pseudo.local + in + if (x > accu) then x + else accu + ) + and lmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + let x = + Positive_int.to_int x.Pseudo.Primitive_non_local.proj + in + if (x > accu) then x + else accu + ) in if (x > accu) then x else accu ) in - if (x > accu) then x - else accu - ) - in - - let () = - Ezfio.set_pseudo_pseudo_klocmax klocmax; - Ezfio.set_pseudo_pseudo_kmax kmax; - Ezfio.set_pseudo_pseudo_lmax lmax; - let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.local ~f:(fun i (y,c) -> - tmp_array_v_k.(i).(j) <- AO_coef.to_float c; - let y, z = - AO_expo.to_float y.Pseudo.Primitive_local.expo, - R_power.to_int y.Pseudo.Primitive_local.r_power + + let kmax = + Array.init (lmax+1) ~f:(fun i-> + List.map pseudo ~f:(fun x -> + List.filter x.Pseudo.non_local ~f:(fun (y,_) -> + (Positive_int.to_int y.Pseudo.Primitive_non_local.proj) = i) + |> List.length ) + |> List.fold ~init:0 ~f:(fun accu x -> + if accu > x then accu else x) + ) + |> Array.fold ~init:0 ~f:(fun accu i -> + if i > accu then i else accu) + in + + + let () = + Ezfio.set_pseudo_pseudo_klocmax klocmax; + Ezfio.set_pseudo_pseudo_kmax kmax; + Ezfio.set_pseudo_pseudo_lmax lmax; + let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 in - tmp_array_dz_k.(i).(j) <- y; - tmp_array_n_k.(i).(j) <- z; - ) - ); - let concat_2d tmp_array = - let data = - Array.map tmp_array ~f:Array.to_list - |> Array.to_list - |> List.concat - in - Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data - in - concat_2d tmp_array_v_k - |> Ezfio.set_pseudo_pseudo_v_k ; - concat_2d tmp_array_dz_k - |> Ezfio.set_pseudo_pseudo_dz_k; - concat_2d tmp_array_n_k - |> Ezfio.set_pseudo_pseudo_n_k; + List.iteri pseudo ~f:(fun j x -> + List.iteri x.Pseudo.local ~f:(fun i (y,c) -> + tmp_array_v_k.(i).(j) <- AO_coef.to_float c; + let y, z = + AO_expo.to_float y.Pseudo.Primitive_local.expo, + R_power.to_int y.Pseudo.Primitive_local.r_power + in + tmp_array_dz_k.(i).(j) <- y; + tmp_array_n_k.(i).(j) <- z; + ) + ); + let concat_2d tmp_array = + let data = + Array.map tmp_array ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data + in + concat_2d tmp_array_v_k + |> Ezfio.set_pseudo_pseudo_v_k ; + concat_2d tmp_array_dz_k + |> Ezfio.set_pseudo_pseudo_dz_k; + concat_2d tmp_array_n_k + |> Ezfio.set_pseudo_pseudo_n_k; - let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 ) - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.non_local ~f:(fun i (y,c) -> - let k, y, z = - Positive_int.to_int y.Pseudo.Primitive_non_local.proj, - AO_expo.to_float y.Pseudo.Primitive_non_local.expo, - R_power.to_int y.Pseudo.Primitive_non_local.r_power + let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 )) in - tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; - tmp_array_dz_kl.(k).(i).(j) <- y; - tmp_array_n_kl.(k).(i).(j) <- z; - ) - ); - let concat_3d tmp_array = - let data = - Array.map tmp_array ~f:(fun x -> - Array.map x ~f:Array.to_list - |> Array.to_list - |> List.concat) - |> Array.to_list + List.iteri pseudo ~f:(fun j x -> + let last_idx = + Array.create ~len:(lmax+1) 0 + in + List.iter x.Pseudo.non_local ~f:(fun (y,c) -> + let k, y, z = + Positive_int.to_int y.Pseudo.Primitive_non_local.proj, + AO_expo.to_float y.Pseudo.Primitive_non_local.expo, + R_power.to_int y.Pseudo.Primitive_non_local.r_power + in + let i = + last_idx.(k) + in + tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; + tmp_array_dz_kl.(k).(i).(j) <- y; + tmp_array_n_kl.(k).(i).(j) <- z; + last_idx.(k) <- i+1; + ) + ); + let concat_3d tmp_array = + let data = + Array.map tmp_array ~f:(fun x -> + Array.map x ~f:Array.to_list + |> Array.to_list + |> List.concat) + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data + in + concat_3d tmp_array_v_kl + |> Ezfio.set_pseudo_pseudo_v_kl ; + concat_3d tmp_array_dz_kl + |> Ezfio.set_pseudo_pseudo_dz_kl ; + concat_3d tmp_array_n_kl + |> Ezfio.set_pseudo_pseudo_n_kl ; + in + + + + + (* Write Basis set *) + let basis = + + let nmax = + Nucl_number.get_max () + in + let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function + | [] -> accu + | e::tail -> + let new_accu = + (e,(Nucl_number.of_int ~max:nmax n))::accu + in + do_work new_accu (n+1) tail + in + let result = do_work [] 1 nuclei + |> List.rev + |> List.map ~f:(fun (x,i) -> + try + let e = + match x.Atom.element with + | Element.X -> Element.H + | e -> e + in + Basis.read_element (basis_channel x.Atom.element) i e + with + | End_of_file -> failwith + ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") + ) |> List.concat + in + (* close all in_channels *) + result in - Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data - in - concat_3d tmp_array_v_kl - |> Ezfio.set_pseudo_pseudo_v_kl ; - concat_3d tmp_array_dz_kl - |> Ezfio.set_pseudo_pseudo_dz_kl ; - concat_3d tmp_array_n_kl - |> Ezfio.set_pseudo_pseudo_n_kl ; - in - - - - - (* Write Basis set *) - let basis = - - let nmax = - Nucl_number.get_max () - in - let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function - | [] -> accu - | e::tail -> - let new_accu = - (e,(Nucl_number.of_int ~max:nmax n))::accu + let long_basis = Long_basis.of_basis basis in + let ao_num = List.length long_basis in + Ezfio.set_ao_basis_ao_num ao_num; + Ezfio.set_ao_basis_ao_basis b; + let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) + and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) + and ao_power= + let l = List.map long_basis ~f:(fun (x,_,_) -> x) in + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) in - do_work new_accu (n+1) tail - in - let result = do_work [] 1 nuclei - |> List.rev - |> List.map ~f:(fun (x,i) -> - try - let e = - match x.Atom.element with - | Element.X -> Element.H - | e -> e - in - Basis.read_element (basis_channel x.Atom.element) i e - with - | End_of_file -> failwith - ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") - ) - |> List.concat - in - (* close all in_channels *) - result - in - let long_basis = Long_basis.of_basis basis in - let ao_num = List.length long_basis in - Ezfio.set_ao_basis_ao_num ao_num; - Ezfio.set_ao_basis_ao_basis b; - let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) - and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) - and ao_power= - let l = List.map long_basis ~f:(fun (x,_,_) -> x) in - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) - in - let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> - if x > s then x - else s) ao_prim_num - in - let gtos = - List.map long_basis ~f:(fun (_,x,_) -> x) - in - - let create_expo_coef ec = - let coefs = - begin match ec with - | `Coefs -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) - | `Expos -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float - prim.Primitive.expo) ) - end + let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> + if x > s then x + else s) ao_prim_num in - let rec get_n n accu = function - | [] -> List.rev accu - | h::tail -> - let y = - begin match List.nth h n with - | Some x -> x - | None -> 0. + let gtos = + List.map long_basis ~f:(fun (_,x,_) -> x) + in + + let create_expo_coef ec = + let coefs = + begin match ec with + | `Coefs -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) + | `Expos -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float + prim.Primitive.expo) ) end - in - get_n n (y::accu) tail + in + let rec get_n n accu = function + | [] -> List.rev accu + | h::tail -> + let y = + begin match List.nth h n with + | Some x -> x + | None -> 0. + end + in + get_n n (y::accu) tail + in + let rec build accu = function + | n when n=ao_prim_num_max -> accu + | n -> build ( accu @ (get_n n [] coefs) ) (n+1) + in + build [] 0 in - let rec build accu = function - | n when n=ao_prim_num_max -> accu - | n -> build ( accu @ (get_n n [] coefs) ) (n+1) - in - build [] 0 - in - let ao_coef = create_expo_coef `Coefs - and ao_expo = create_expo_coef `Expos + let ao_coef = create_expo_coef `Coefs + and ao_expo = create_expo_coef `Expos + in + let () = + Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; + Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; + Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; + Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; + Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; + Ezfio.set_ao_basis_ao_cartesian(cart); + in + match Input.Ao_basis.read () with + | None -> failwith "Error in basis" + | Some x -> Input.Ao_basis.write x in let () = - Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; - Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; - Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; - Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; - Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; - Ezfio.set_ao_basis_ao_cartesian(cart); - in - match Input.Ao_basis.read () with - | None -> failwith "Error in basis" - | Some x -> Input.Ao_basis.write x + try write_file () with + | ex -> + begin + begin + match Sys.is_directory ezfio_file with + | `Yes -> rmdir ezfio_file + | _ -> () + end; + raise ex; + end + in ()