From 4beca8f42dce9add8088c0f0ddae49c035e3286e Mon Sep 17 00:00:00 2001 From: Manu Date: Fri, 24 Oct 2014 19:28:54 +0200 Subject: [PATCH] Starting work on DDCI --- ocaml/qp_set_ddci.ml | 303 ++++++++++++++++++++++++++++++ src/Bitmask/bitmasks.ezfio_config | 2 +- src/Bitmask/bitmasks.irp.f | 4 - 3 files changed, 304 insertions(+), 5 deletions(-) create mode 100644 ocaml/qp_set_ddci.ml diff --git a/ocaml/qp_set_ddci.ml b/ocaml/qp_set_ddci.ml new file mode 100644 index 00000000..aa6aedd4 --- /dev/null +++ b/ocaml/qp_set_ddci.ml @@ -0,0 +1,303 @@ +open Qputils;; +open Qptypes;; +open Core.Std;; + +(* + * Command-line arguments + * ---------------------- + *) + +let build_mask from upto n_int = + let from = MO_number.to_int from + and upto = MO_number.to_int upto + and n_int = N_int_number.to_int n_int + in + let rec build_mask bit = function + | 0 -> [] + | i -> + if ( i = upto ) then + Bit.One::(build_mask Bit.One (i-1)) + else if ( i = from ) then + Bit.One::(build_mask Bit.Zero (i-1)) + else + bit::(build_mask bit (i-1)) + in + let starting_bit = + if ( (upto >= n_int*64) || (upto < 0) ) then Bit.One + else Bit.Zero + in + build_mask starting_bit (n_int*64) + |> List.rev +;; + + +let failure s = raise (Failure s) +;; + +type t = + | Core + | Inactive + | Active + | Virtual + | Deleted + | None +;; + +let t_to_string = function + | Core -> "core" + | Inactive -> "inactive" + | Active -> "active" + | Virtual -> "virtual" + | Deleted -> "deleted" + | None -> assert false +;; + +let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename = + + Ezfio.set_file ezfio_filename ; + if not (Ezfio.has_mo_basis_mo_tot_num ()) then + failure "mo_basis/mo_tot_num not found" ; + + let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in + let n_int = + try N_int_number.of_int (Ezfio.get_determinants_n_int ()) + with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num + in + + + let mo_class = Array.init mo_tot_num ~f:(fun i -> None) in + + (* Check input data *) + let apply_class l = + let rec apply_class t = function + | [] -> () + | k::tail -> let i = MO_number.to_int k in + begin + match mo_class.(i-1) with + | None -> mo_class.(i-1) <- t ; + apply_class t tail; + | x -> failure + (Printf.sprintf "Orbital %d is defined both in the %s and %s spaces" + i (t_to_string x) (t_to_string t)) + end + in + match l with + | MO_class.Core x -> apply_class Core x + | MO_class.Inactive x -> apply_class Inactive x + | MO_class.Active x -> apply_class Active x + | MO_class.Virtual x -> apply_class Virtual x + | MO_class.Deleted x -> apply_class Deleted x + in + + let core_input = core in + let core = MO_class.create_core core in + let inact = MO_class.create_inactive inact in + let act = MO_class.create_active act in + let virt = MO_class.create_virtual virt in + let del = MO_class.create_deleted del in + + apply_class core ; + apply_class inact ; + apply_class act ; + apply_class virt ; + apply_class del ; + + for i=1 to (Array.length mo_class) + do + if (mo_class.(i-1) = None) then + failure (Printf.sprintf "Orbital %d is not specified (mo_tot_num = %d)" i mo_tot_num) + done; + + + (* Debug output *) + MO_class.to_string core |> print_endline ; + MO_class.to_string inact |> print_endline ; + MO_class.to_string act |> print_endline ; + MO_class.to_string virt |> print_endline ; + MO_class.to_string del |> print_endline ; + + (* Create masks *) + let ia = Excitation.create_single inact act + and aa = Excitation.create_single act act + and av = Excitation.create_single act virt + and iv = Excitation.create_single inact virt + in + let single_excitations = [ ia ; aa ; av ; iv] + |> List.map ~f:Excitation.(fun x -> + match x with + | Single (x,y) -> + ( MO_class.to_bitlist n_int (Hole.to_mo_class x), + MO_class.to_bitlist n_int (Particle.to_mo_class y) ) + | Double _ -> assert false + ) + + and double_excitations = [ + Excitation.double_of_singles ia ia ; + Excitation.double_of_singles ia aa ; + Excitation.double_of_singles ia iv ; + Excitation.double_of_singles ia av ; + + Excitation.double_of_singles aa aa ; + Excitation.double_of_singles aa iv ; + Excitation.double_of_singles aa av ; + + Excitation.double_of_singles iv aa ; + Excitation.double_of_singles iv av ; + + Excitation.double_of_singles iv iv ] + + |> List.map ~f:Excitation.(fun x -> + match x with + | Single _ -> assert false + | Double (x,y,z,t) -> + ( MO_class.to_bitlist n_int (Hole.to_mo_class x), + MO_class.to_bitlist n_int (Particle.to_mo_class y) , + MO_class.to_bitlist n_int (Hole.to_mo_class z), + MO_class.to_bitlist n_int (Particle.to_mo_class t) ) + ) + in + + let extract_hole (h,_) = h + and extract_particle (_,p) = p + and extract_hole1 (h,_,_,_) = h + and extract_particle1 (_,p,_,_) = p + and extract_hole2 (_,_,h,_) = h + and extract_particle2 (_,_,_,p) = p + in + let result_ref = + let core = MO_class.create_inactive core_input in + let cv = Excitation.create_single core virt in + let cv = match cv with + | Single (x,y) -> + ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x), + MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) + | Double _ -> assert false + in + let iv = match iv with + | Single (x,y) -> + ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x), + MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) + | Double _ -> assert false + in + [ Bitlist.or_operator (extract_hole iv) (extract_hole cv); + extract_particle iv ] + in + + let result_gen = [ + List.map ~f:extract_hole single_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + List.map ~f:extract_particle single_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + List.map ~f:extract_hole1 double_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + List.map ~f:extract_particle1 double_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + List.map ~f:extract_hole2 double_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + List.map ~f:extract_particle2 double_excitations + |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + + ] + in + + (* Print bitmasks *) + print_endline "Reference Bitmasks:"; + List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_ref; + + print_endline "Generators Bitmasks:"; + List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_gen; + + (* Transform to int64 *) + let result_gen = List.map ~f:(fun x -> + let y = Bitlist.to_int64_list x in y@y ) + result_gen + |> List.concat + in + let result_ref = List.map ~f:(fun x -> + let y = Bitlist.to_int64_list x in y@y ) + result_ref + |> List.concat + in + + (* Write generators masks *) + Ezfio.set_bitmasks_n_int (N_int_number.to_int n_int); + Ezfio.set_bitmasks_bit_kind 8; + Ezfio.set_bitmasks_n_mask_gen 1; + Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 6; 1|]) ~data:result_gen + |> Ezfio.set_bitmasks_generators ; + + (* Write reference masks *) + Ezfio.set_bitmasks_n_mask_ref 1; + Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 2; 1|]) ~data:result_ref + |> Ezfio.set_bitmasks_reference ; + + +;; + +let ezfio_file = + let failure filename = + eprintf "'%s' is not an EZFIO file.\n%!" filename; + exit 1 + in + Command.Spec.Arg_type.create + (fun filename -> + match Sys.is_directory filename with + | `Yes -> + begin + match Sys.is_file (filename ^ "/.version") with + | `Yes -> filename + | _ -> failure filename + end + | _ -> failure filename + ) +;; + +let default range = + let failure filename = + eprintf "'%s' is not a regular file.\n%!" filename; + exit 1 + in + Command.Spec.Arg_type.create + (fun filename -> + match Sys.is_directory filename with + | `Yes -> + begin + match Sys.is_file (filename^"/.version") with + | `Yes -> filename + | _ -> failure filename + end + | _ -> failure filename + ) +;; + +let spec = + let open Command.Spec in + empty + +> flag "core" (optional string) ~doc:"range Range of core orbitals" + +> flag "inact" (optional string) ~doc:"range Range of inactive orbitals" + +> flag "act" (optional string) ~doc:"range Range of active orbitals" + +> flag "virt" (optional string) ~doc:"range Range of virtual orbitals" + +> flag "del" (optional string) ~doc:"range Range of deleted orbitals" + +> anon ("ezfio_filename" %: ezfio_file) +;; + +let command = + Command.basic + ~summary: "Quantum Package command" + ~readme:(fun () -> + "Set the orbital classes in an EZFIO directory + The range of MOs has the form : \"[36-53,72-107,126-131]\" + ") + spec + (fun core inact act virt del ezfio_filename () -> run ?core ?inact ?act ?virt ?del ezfio_filename ) +;; + +let () = + Command.run command + + diff --git a/src/Bitmask/bitmasks.ezfio_config b/src/Bitmask/bitmasks.ezfio_config index 79f1ae6c..988510fb 100644 --- a/src/Bitmask/bitmasks.ezfio_config +++ b/src/Bitmask/bitmasks.ezfio_config @@ -4,5 +4,5 @@ bitmasks N_mask_gen integer generators integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen) N_mask_ref integer - reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_ref) + reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,2,bitmasks_N_mask_ref) diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 353d6a3b..dbf3304c 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -175,10 +175,6 @@ BEGIN_PROVIDER [ integer(bit_kind), reference_bitmask, (N_int,2,2,N_reference_bi else reference_bitmask(:,:,s_hole ,1) = HF_bitmask reference_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - reference_bitmask(:,:,d_hole1,1) = HF_bitmask - reference_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - reference_bitmask(:,:,d_hole2,1) = HF_bitmask - reference_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) endif END_PROVIDER