diff --git a/TODO b/TODO index 046510ed..abdb618f 100644 --- a/TODO +++ b/TODO @@ -2,16 +2,8 @@ * Faire que le slave de Hartree-fock est le calcul des integrales AO en parallele -# Web/doc - -* Creer une page web pas trop degueu et la mettre ici : http://lcpq.github.io/quantum_package - -* Creer une page avec la liste de tous les exectuables - - # Exterieur -* Molden format : http://cheminf.cmbi.ru.nl/molden/molden_format.html : read+write. Thomas est dessus * Un module pour lire les integrales Moleculaires depuis un FCIDUMP * Un module pour lire des integrales Atomiques (voir module de Mimi pour lire les AO Slater) * Format Fchk (gaussian) @@ -24,51 +16,22 @@ # User doc: - * Videos: - +) RHF - * Renvoyer a la doc des modules : c'est pour les programmeurs au depart! * Mettre le mp2 comme exercice - * Interfaces : molden/fcidump - * Natural orbitals - * Parameters for Hartree-Fock - * Parameters for Davidson - * Running in parallel # Programmers doc: * Example : Simple Hartree-Fock program from scratch * Examples : subroutine example_module +# enleverle psi_det_size for all complicated stuffs with dimension of psi_coef + # Config file for Cray -# EZFIO sans fork - -Refaire les benchmarks - -# Documentation de qpsh - # Documentation de /etc -# Toto -Re-design de qp command - -Doc: plugins et qp_plugins - Ajouter les symetries dans devel -<<<<<<< HEAD -Compiler ezfio avec openmp - -# Parallelize i_H_psi -======= - -# Parallelize i_H_psi -<<<<<<< HEAD -======= - - ->>>>>>> minor_modifs IMPORTANT: Davidson Diagonalization diff --git a/configure b/configure index 21df1b71..aa27ffa4 100755 --- a/configure +++ b/configure @@ -3,18 +3,23 @@ # Quantum Package configuration script # +unset CC +unset CXX + TEMP=$(getopt -o c:i:h -l config:,install:,help -n $0 -- "$@") || exit 1 eval set -- "$TEMP" export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )" echo "QP_ROOT="$QP_ROOT +unset CC +unset CCXX # When updating version, update also etc files BATS_URL="https://github.com/bats-core/bats-core/archive/v1.1.0.tar.gz" BUBBLE_URL="https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz" DOCOPT_URL="https://github.com/docopt/docopt/archive/0.6.2.tar.gz" -EZFIO_URL="https://gitlab.com/scemama/EZFIO/-/archive/v1.4.0/EZFIO-v1.4.0.tar.gz" +EZFIO_URL="https://gitlab.com/scemama/EZFIO/-/archive/v1.6.1/EZFIO-v1.6.1.tar.gz" F77ZMQ_URL="https://github.com/scemama/f77_zmq/archive/v4.2.5.tar.gz" GMP_URL="ftp://ftp.gnu.org/gnu/gmp/gmp-6.1.2.tar.bz2" IRPF90_URL="https://gitlab.com/scemama/irpf90/-/archive/v1.7.6/irpf90-v1.7.6.tar.gz" diff --git a/data/basis/aug-cc-pvtz b/data/basis/aug-cc-pvtz index b9d1788f..5a5fd369 100644 --- a/data/basis/aug-cc-pvtz +++ b/data/basis/aug-cc-pvtz @@ -92,52 +92,58 @@ F 1 1 0.0816000 1.0000000 BERYLLIUM -S 9 - 1 6863.0000000 0.0002360 - 2 1030.0000000 0.0018260 - 3 234.7000000 0.0094520 - 4 66.5600000 0.0379570 - 5 21.6900000 0.1199650 - 6 7.7340000 0.2821620 - 7 2.9160000 0.4274040 - 8 1.1300000 0.2662780 - 9 0.1101000 -0.0072750 -S 9 - 1 6863.0000000 -0.0000430 - 2 1030.0000000 -0.0003330 - 3 234.7000000 -0.0017360 - 4 66.5600000 -0.0070120 - 5 21.6900000 -0.0231260 - 6 7.7340000 -0.0581380 - 7 2.9160000 -0.1145560 - 8 1.1300000 -0.1359080 - 9 0.1101000 0.5774410 +S 11 +1 6.863000E+03 2.360000E-04 +2 1.030000E+03 1.826000E-03 +3 2.347000E+02 9.452000E-03 +4 6.656000E+01 3.795700E-02 +5 2.169000E+01 1.199650E-01 +6 7.734000E+00 2.821620E-01 +7 2.916000E+00 4.274040E-01 +8 1.130000E+00 2.662780E-01 +9 2.577000E-01 1.819300E-02 +10 1.101000E-01 -7.275000E-03 +11 4.409000E-02 1.903000E-03 +S 11 +1 6.863000E+03 -4.300000E-05 +2 1.030000E+03 -3.330000E-04 +3 2.347000E+02 -1.736000E-03 +4 6.656000E+01 -7.012000E-03 +5 2.169000E+01 -2.312600E-02 +6 7.734000E+00 -5.813800E-02 +7 2.916000E+00 -1.145560E-01 +8 1.130000E+00 -1.359080E-01 +9 2.577000E-01 2.280260E-01 +10 1.101000E-01 5.774410E-01 +11 4.409000E-02 3.178730E-01 S 1 - 1 0.2577000 1.0000000 +1 2.577000E-01 1.000000E+00 S 1 - 1 0.0440900 1.0000000 +1 4.409000E-02 1.000000E+00 S 1 - 1 0.0150300 1.0000000 -P 3 - 1 7.4360000 0.0107360 - 2 1.5770000 0.0628540 - 3 0.4352000 0.2481800 +1 1.470000E-02 1.000000E+00 +P 5 +1 7.436000E+00 1.073600E-02 +2 1.577000E+00 6.285400E-02 +3 4.352000E-01 2.481800E-01 +4 1.438000E-01 5.236990E-01 +5 4.994000E-02 3.534250E-01 P 1 - 1 0.1438000 1.0000000 +1 1.438000E-01 1.000000E+00 P 1 - 1 0.0499400 1.0000000 +1 4.994000E-02 1.000000E+00 P 1 - 1 0.0070600 1.0000000 +1 9.300000E-03 1.000000E+00 D 1 - 1 0.3480000 1.0000000 +1 3.493000E-01 1.000000E+00 D 1 - 1 0.1803000 1.0000000 +1 1.724000E-01 1.000000E+00 D 1 - 1 0.0654000 1.0000000 +1 5.880000E-02 1.000000E+00 F 1 - 1 0.3250000 1.0000000 +1 3.423000E-01 1.0000000 F 1 - 1 0.1533000 1.0000000 +1 1.188000E-01 1.000000E+00 BORON S 8 diff --git a/ocaml/Input_bitmasks.ml b/ocaml/Input_bitmasks.ml index 944a80ff..921b34da 100644 --- a/ocaml/Input_bitmasks.ml +++ b/ocaml/Input_bitmasks.ml @@ -6,10 +6,6 @@ module Bitmasks : sig type t = { n_int : N_int_number.t; bit_kind : Bit_kind.t; - n_mask_gen : Bitmask_number.t; - generators : int64 array; - n_mask_cas : Bitmask_number.t; - cas : int64 array; } [@@deriving sexp] ;; val read : unit -> t option @@ -18,12 +14,7 @@ end = struct type t = { n_int : N_int_number.t; bit_kind : Bit_kind.t; - n_mask_gen : Bitmask_number.t; - generators : int64 array; - n_mask_cas : Bitmask_number.t; - cas : int64 array; } [@@deriving sexp] - ;; let get_default = Qpackage.get_ezfio_default "bitmasks";; @@ -36,7 +27,6 @@ end = struct ; Ezfio.get_bitmasks_n_int () |> N_int_number.of_int - ;; let read_bit_kind () = if not (Ezfio.has_bitmasks_bit_kind ()) then @@ -46,89 +36,12 @@ end = struct ; Ezfio.get_bitmasks_bit_kind () |> Bit_kind.of_int - ;; - - let read_n_mask_gen () = - if not (Ezfio.has_bitmasks_n_mask_gen ()) then - Ezfio.set_bitmasks_n_mask_gen 1 - ; - Ezfio.get_bitmasks_n_mask_gen () - |> Bitmask_number.of_int - ;; - - - let full_mask n_int = - let range = "[1-"^ - (string_of_int (Ezfio.get_mo_basis_mo_num ()))^"]" - in - MO_class.create_active range - |> MO_class.to_bitlist n_int - ;; - - let read_generators () = - if not (Ezfio.has_bitmasks_generators ()) then - begin - let n_int = - read_n_int () - in - let act = - full_mask n_int - in - let result = [ act ; act ; act ; act ; act ; act ] - |> List.map (fun x -> - let y = Bitlist.to_int64_list x in y@y ) - |> List.concat - in - let generators = Ezfio.ezfio_array_of_list ~rank:4 - ~dim:([| (N_int_number.to_int n_int) ; 2; 6; 1|]) ~data:result - in - Ezfio.set_bitmasks_generators generators - end; - Ezfio.get_bitmasks_generators () - |> Ezfio.flattened_ezfio - ;; - - let read_n_mask_cas () = - if not (Ezfio.has_bitmasks_n_mask_cas ()) then - Ezfio.set_bitmasks_n_mask_cas 1 - ; - Ezfio.get_bitmasks_n_mask_cas () - |> Bitmask_number.of_int - ;; - - - let read_cas () = - if not (Ezfio.has_bitmasks_cas ()) then - begin - let n_int = - read_n_int () - in - let act = - full_mask n_int - in - let result = [ act ; act ] - |> List.map (fun x -> - let y = Bitlist.to_int64_list x in y@y ) - |> List.concat - in - let cas = Ezfio.ezfio_array_of_list ~rank:3 - ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result - in - Ezfio.set_bitmasks_cas cas - end; - Ezfio.get_bitmasks_cas () - |> Ezfio.flattened_ezfio - ;; let read () = if (Ezfio.has_mo_basis_mo_num ()) then Some { n_int = read_n_int (); bit_kind = read_bit_kind (); - n_mask_gen = read_n_mask_gen (); - generators = read_generators (); - n_mask_cas = read_n_mask_cas (); - cas = read_cas (); } else None @@ -138,21 +51,9 @@ end = struct Printf.sprintf " n_int = %s bit_kind = %s -n_mask_gen = %s -generators = %s -n_mask_cas = %s -cas = %s " (N_int_number.to_string b.n_int) (Bit_kind.to_string b.bit_kind) - (Bitmask_number.to_string b.n_mask_gen) - (Array.to_list b.generators - |> List.map (fun x-> Int64.to_string x) - |> String.concat ", ") - (Bitmask_number.to_string b.n_mask_cas) - (Array.to_list b.cas - |> List.map (fun x-> Int64.to_string x) - |> String.concat ", ") end diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index 6c449c1b..9c316f8c 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -15,7 +15,7 @@ module Determinants_by_hand : sig state_average_weight : Positive_float.t array; } [@@deriving sexp] val read : ?full:bool -> unit -> t option - val write : t -> unit + val write : ?force:bool -> t -> unit val to_string : t -> string val to_rst : t -> Rst_string.t val of_rst : Rst_string.t -> t option @@ -318,22 +318,23 @@ end = struct None ;; - let write { n_int ; - bit_kind ; - n_det ; - n_det_qp_edit ; - expected_s2 ; - psi_coef ; - psi_det ; - n_states ; - state_average_weight ; - } = + let write ?(force=false) + { n_int ; + bit_kind ; + n_det ; + n_det_qp_edit ; + expected_s2 ; + psi_coef ; + psi_det ; + n_states ; + state_average_weight ; + } = write_n_int n_int ; write_bit_kind bit_kind; write_n_det n_det; write_n_states n_states; write_expected_s2 expected_s2; - if n_det <= n_det_qp_edit then + if force || (n_det <= n_det_qp_edit) then begin write_n_det_qp_edit n_det; write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ; @@ -596,7 +597,7 @@ psi_det = %s let new_det = { det with n_det = (Det_number.of_int n_det_new) } in - write new_det + write ~force:true new_det ;; let extract_state istate = @@ -628,7 +629,7 @@ psi_det = %s let new_det = { det with n_states = (States_number.of_int 1) } in - write new_det + write ~force:true new_det ;; let extract_states range = @@ -665,6 +666,7 @@ psi_det = %s det.psi_coef.(!state_shift+i) <- det.psi_coef.(i+ishift) done + ; Printf.printf "OK\n%!" ; end; state_shift := !state_shift + n_det ) sorted_list @@ -672,7 +674,7 @@ psi_det = %s let new_det = { det with n_states = (States_number.of_int @@ List.length sorted_list) } in - write new_det + write ~force:true new_det ;; end diff --git a/ocaml/Input_nuclei_by_hand.ml b/ocaml/Input_nuclei_by_hand.ml index 520b4f05..f195a2de 100644 --- a/ocaml/Input_nuclei_by_hand.ml +++ b/ocaml/Input_nuclei_by_hand.ml @@ -175,7 +175,7 @@ nucl_coord = %s nucl_num ) :: ( List.init nucl_num (fun i-> - Printf.sprintf " %-3s %d %s" + Printf.sprintf " %-3s %3d %s" (b.nucl_label.(i) |> Element.to_string) (b.nucl_charge.(i) |> Charge.to_int ) (b.nucl_coord.(i) |> Point3d.to_string ~units:Units.Angstrom) ) diff --git a/ocaml/Makefile b/ocaml/Makefile index 6ff91273..978f7e87 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -80,7 +80,7 @@ git: ./create_git_sha1.sh ${QP_EZFIO}/Ocaml/ezfio.ml: - $(NINJA) -C ${QP_EZFIO} + $(NINJA) -C ${QP_ROOT}/config ${QP_ROOT}/lib/libezfio_irp.a qp_edit.ml: ../scripts/ezfio_interface/qp_edit_template diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index 942e2cc2..7ab861e2 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/ocaml/qp_set_mo_class.ml @@ -106,96 +106,6 @@ let set ~core ~inact ~act ~virt ~del = 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 - in - let single_excitations = [ ia ; aa ; av ] - |> List.map (fun z -> - let open Excitation in - match z 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 av ; - Excitation.double_of_singles aa aa ; - Excitation.double_of_singles aa av ; - Excitation.double_of_singles av av ] - |> List.map (fun x -> - let open Excitation in - 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 init = Bitlist.zero n_int in - let result = [ - List.map extract_hole single_excitations - |> List.fold_left Bitlist.or_operator init; - List.map extract_particle single_excitations - |> List.fold_left Bitlist.or_operator init; - List.map extract_hole1 double_excitations - |> List.fold_left Bitlist.or_operator init; - List.map extract_particle1 double_excitations - |> List.fold_left Bitlist.or_operator init; - List.map extract_hole2 double_excitations - |> List.fold_left Bitlist.or_operator init; - List.map extract_particle2 double_excitations - |> List.fold_left Bitlist.or_operator init; - ] - in - - (* Debug masks in output - List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result; - *) - - (* Write masks *) - let result = - List.map (fun x -> - let y = Bitlist.to_int64_list x in y@y ) - result - |> List.concat - in - - 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 - |> Ezfio.set_bitmasks_generators ; - - let result = - let open Excitation in - match aa with - | Double _ -> assert false - | Single (x,y) -> - Bitlist.to_int64_list - ( MO_class.to_bitlist n_int ( Hole.to_mo_class x) ) @ - Bitlist.to_int64_list - ( MO_class.to_bitlist n_int (Particle.to_mo_class y) ) - in - Ezfio.set_bitmasks_n_mask_cas 1; - Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result - |> Ezfio.set_bitmasks_cas; - let data = Array.to_list mo_class |> List.map (fun x -> match x with diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index a63a19cc..2c54a218 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -78,9 +78,6 @@ let input_data = " | _ -> raise (Invalid_argument \"Bit_kind should be (1|2|4|8).\") end; -* Bitmask_number : int - assert (x > 0) ; - * MO_coef : float * MO_occ : float diff --git a/scripts/compilation/qp_create_ninja b/scripts/compilation/qp_create_ninja index 7a148773..8381d1a2 100755 --- a/scripts/compilation/qp_create_ninja +++ b/scripts/compilation/qp_create_ninja @@ -839,21 +839,6 @@ if __name__ == "__main__": l_module = d_binaries.keys() - # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # - # C h e c k _ c o h e r e n c y # - # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # - - for module in dict_root_path.values(): - - if module not in d_binaries: - l_msg = ["{0} is a root module but does not contain a main file.", - "- Create it in {0}", - "- Or delete {0} `qp_module uninstall {0}`", - "- Or install a module that needs {0} with a main "] - - print "\n".join(l_msg).format(module.rel) - sys.exit(1) - # ~#~#~#~#~#~#~#~#~#~#~#~ # # G l o b a l _ b u i l d # # ~#~#~#~#~#~#~#~#~#~#~#~ # diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index d7c3fd32..4218456d 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -120,7 +120,7 @@ let set str s = match s with {write} | Electrons -> write Electrons.(of_rst, write) s - | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s + | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write ~force:false) s | Nuclei_by_hand -> write Nuclei_by_hand.(of_rst, write) s | Ao_basis -> () (* TODO *) | Mo_basis -> () (* TODO *) diff --git a/src/bitmask/bitmask_cas_routines.irp.f b/src/bitmask/bitmask_cas_routines.irp.f index c0c8cd11..4c3faebe 100644 --- a/src/bitmask/bitmask_cas_routines.irp.f +++ b/src/bitmask/bitmask_cas_routines.irp.f @@ -3,28 +3,28 @@ integer function number_of_holes(key_in) BEGIN_DOC ! Function that returns the number of holes in the inact space ! -! popcnt( -! xor( -! iand( -! reunion_of_core_inact_bitmask(1,1), -! xor( -! key_in(1,1), -! iand( -! key_in(1,1), -! cas_bitmask(1,1,1)) -! ) -! ), -! reunion_of_core_inact_bitmask(1,1)) ) -! -! (key_in && cas_bitmask) -! +---------------------+ -! electrons in cas xor key_in -! +---------------------------------+ -! electrons outside of cas && reunion_of_core_inact_bitmask -! +------------------------------------------------------------------+ -! electrons in the core/inact space xor reunion_of_core_inact_bitmask -! +---------------------------------------------------------------------------------+ -! holes + ! popcnt( + ! xor( + ! iand( + ! reunion_of_core_inact_bitmask(1,1), + ! xor( + ! key_in(1,1), + ! iand( + ! key_in(1,1), + ! act_bitmask(1,1)) + ! ) + ! ), + ! reunion_of_core_inact_bitmask(1,1)) ) + ! + ! (key_in && act_bitmask) + ! +---------------------+ + ! electrons in cas xor key_in + ! +---------------------------------+ + ! electrons outside of cas && reunion_of_core_inact_bitmask + ! +------------------------------------------------------------------+ + ! electrons in the core/inact space xor reunion_of_core_inact_bitmask + ! +---------------------------------------------------------------------------------+ + ! holes END_DOC implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) @@ -33,74 +33,32 @@ integer function number_of_holes(key_in) if(N_int == 1)then number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) else if(N_int == 2)then number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) ) else if(N_int == 3)then number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1)))), reunion_of_core_inact_bitmask(3,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2)))), reunion_of_core_inact_bitmask(3,2)) ) else if(N_int == 4)then number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) - else if(N_int == 5)then - number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) ) - else if(N_int == 6)then - number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) ) - else if(N_int == 7)then - number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1)))), reunion_of_core_inact_bitmask(3,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2)))), reunion_of_core_inact_bitmask(3,2)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),act_bitmask(4,1)))), reunion_of_core_inact_bitmask(4,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),act_bitmask(4,2)))), reunion_of_core_inact_bitmask(4,2)) ) else do i = 1, N_int number_of_holes = number_of_holes & @@ -111,11 +69,11 @@ integer function number_of_holes(key_in) xor( & key_in(i,1), & ! MOs of key_in not in the CAS iand( & ! MOs of key_in in the CAS - key_in(i,1), cas_bitmask(i,1,1) & + key_in(i,1), act_bitmask(i,1) & ) & ) & ), reunion_of_core_inact_bitmask(i,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1)))), reunion_of_core_inact_bitmask(i,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),act_bitmask(i,2)))), reunion_of_core_inact_bitmask(i,2)) ) enddo endif end @@ -131,97 +89,37 @@ integer function number_of_particles(key_in) number_of_particles= 0 if(N_int == 1)then number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) )) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) )) else if(N_int == 2)then number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) ) ) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) ) ) else if(N_int == 3)then number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) )) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) )) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) )) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) )) & + + popcnt( iand( xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1))), virt_bitmask(3,1) )) & + + popcnt( iand( xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2))), virt_bitmask(3,2) )) else if(N_int == 4)then number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) - else if(N_int == 5)then - number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) - else if(N_int == 6)then - number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) - else if(N_int == 7)then - number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) & - + popcnt( iand( iand( xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1))), virt_bitmask(7,1) ), virt_bitmask(7,1)) ) & - + popcnt( iand( iand( xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1))), virt_bitmask(7,2) ), virt_bitmask(7,2)) ) - else if(N_int == 8)then - number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) & - + popcnt( iand( iand( xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1))), virt_bitmask(7,1) ), virt_bitmask(7,1)) ) & - + popcnt( iand( iand( xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1))), virt_bitmask(7,2) ), virt_bitmask(7,2)) ) & - + popcnt( iand( iand( xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1))), virt_bitmask(8,1) ), virt_bitmask(8,1)) ) & - + popcnt( iand( iand( xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1))), virt_bitmask(8,2) ), virt_bitmask(8,2)) ) + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) ) ) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) ) ) & + + popcnt( iand( xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1))), virt_bitmask(3,1) ) ) & + + popcnt( iand( xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2))), virt_bitmask(3,2) ) ) & + + popcnt( iand( xor(key_in(4,1),iand(key_in(4,1),act_bitmask(4,1))), virt_bitmask(4,1) ) ) & + + popcnt( iand( xor(key_in(4,2),iand(key_in(4,2),act_bitmask(4,2))), virt_bitmask(4,2) ) ) else do i = 1, N_int - number_of_particles= number_of_particles & - + popcnt( iand( iand( xor(key_in(i,1),iand(key_in(i,1),cas_bitmask(i,1,1))), virt_bitmask(i,1) ), virt_bitmask(i,1)) ) & - + popcnt( iand( iand( xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1))), virt_bitmask(i,2) ), virt_bitmask(i,2)) ) + number_of_particles= number_of_particles & + + popcnt( iand( xor(key_in(i,1),iand(key_in(i,1),act_bitmask(i,1))), virt_bitmask(i,1) )) & + + popcnt( iand( xor(key_in(i,2),iand(key_in(i,2),act_bitmask(i,2))), virt_bitmask(i,2) )) enddo endif end @@ -230,7 +128,7 @@ logical function is_a_two_holes_two_particles(key_in) BEGIN_DOC ! logical function that returns True if the determinant 'key_in' ! belongs to the 2h-2p excitation class of the DDCI space - ! this is calculated using the CAS_bitmask that defines the active + ! this is calculated using the act_bitmask that defines the active ! orbital space, the inact_bitmasl that defines the inactive oribital space ! and the virt_bitmask that defines the virtual orbital space END_DOC @@ -246,174 +144,62 @@ logical function is_a_two_holes_two_particles(key_in) i_diff = 0 if(N_int == 1)then i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) & + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) else if(N_int == 2)then i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) & + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) ) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) )) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) )) else if(N_int == 3)then i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) & + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) ) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) ) ) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1)))), reunion_of_core_inact_bitmask(3,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2)))), reunion_of_core_inact_bitmask(3,2)) ) & + + popcnt( iand( xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1))), virt_bitmask(3,1) ) ) & + + popcnt( iand( xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2))), virt_bitmask(3,2) ) ) else if(N_int == 4)then i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) - else if(N_int == 5)then - i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) - else if(N_int == 6)then - i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) - else if(N_int == 7)then - i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) ) & - + popcnt( iand( iand( xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1))), virt_bitmask(7,1) ), virt_bitmask(7,1)) ) & - + popcnt( iand( iand( xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1))), virt_bitmask(7,2) ), virt_bitmask(7,2)) ) - else if(N_int == 8)then - i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) ) & - + popcnt( iand( iand( xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1))), virt_bitmask(2,1) ), virt_bitmask(2,1)) ) & - + popcnt( iand( iand( xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1))), virt_bitmask(2,2) ), virt_bitmask(2,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) ) & - + popcnt( iand( iand( xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1))), virt_bitmask(3,1) ), virt_bitmask(3,1)) ) & - + popcnt( iand( iand( xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1))), virt_bitmask(3,2) ), virt_bitmask(3,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) ) & - + popcnt( iand( iand( xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1))), virt_bitmask(4,1) ), virt_bitmask(4,1)) ) & - + popcnt( iand( iand( xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1))), virt_bitmask(4,2) ), virt_bitmask(4,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) ) & - + popcnt( iand( iand( xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1))), virt_bitmask(5,1) ), virt_bitmask(5,1)) ) & - + popcnt( iand( iand( xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1))), virt_bitmask(5,2) ), virt_bitmask(5,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) ) & - + popcnt( iand( iand( xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1))), virt_bitmask(6,1) ), virt_bitmask(6,1)) ) & - + popcnt( iand( iand( xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1))), virt_bitmask(6,2) ), virt_bitmask(6,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) ) & - + popcnt( iand( iand( xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1))), virt_bitmask(7,1) ), virt_bitmask(7,1)) ) & - + popcnt( iand( iand( xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1))), virt_bitmask(7,2) ), virt_bitmask(7,2)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,1), xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1)))), reunion_of_core_inact_bitmask(8,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,2), xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1)))), reunion_of_core_inact_bitmask(8,2)) ) & - + popcnt( iand( iand( xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1))), virt_bitmask(8,1) ), virt_bitmask(8,1)) ) & - + popcnt( iand( iand( xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1))), virt_bitmask(8,2) ), virt_bitmask(8,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) & + + popcnt( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ) ) & + + popcnt( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1)))), reunion_of_core_inact_bitmask(2,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2)))), reunion_of_core_inact_bitmask(2,2)) ) & + + popcnt( iand( xor(key_in(2,1),iand(key_in(2,1),act_bitmask(2,1))), virt_bitmask(2,1) ) ) & + + popcnt( iand( xor(key_in(2,2),iand(key_in(2,2),act_bitmask(2,2))), virt_bitmask(2,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1)))), reunion_of_core_inact_bitmask(3,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),act_bitmask(3,2)))), reunion_of_core_inact_bitmask(3,2)) ) & + + popcnt( iand( xor(key_in(3,1),iand(key_in(3,1),act_bitmask(3,1))), virt_bitmask(3,1) ) ) & + + popcnt( iand( xor(key_in(4,2),iand(key_in(3,2),act_bitmask(3,2))), virt_bitmask(3,2) ) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),act_bitmask(4,1)))), reunion_of_core_inact_bitmask(4,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),act_bitmask(4,2)))), reunion_of_core_inact_bitmask(4,2)) ) & + + popcnt( iand( xor(key_in(4,1),iand(key_in(4,1),act_bitmask(4,1))), virt_bitmask(4,1) ) ) & + + popcnt( iand( xor(key_in(4,2),iand(key_in(4,2),act_bitmask(4,2))), virt_bitmask(4,2) ) ) else do i = 1, N_int i_diff = i_diff & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,1),iand(key_in(i,1),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1)))), reunion_of_core_inact_bitmask(i,2)) ) & - + popcnt( iand( iand( xor(key_in(i,1),iand(key_in(i,1),cas_bitmask(i,1,1))), virt_bitmask(i,1) ), virt_bitmask(i,1)) ) & - + popcnt( iand( iand( xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1))), virt_bitmask(i,2) ), virt_bitmask(i,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,1),iand(key_in(i,1),act_bitmask(i,1)))), reunion_of_core_inact_bitmask(i,1)) ) & + + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),act_bitmask(i,2)))), reunion_of_core_inact_bitmask(i,2)) ) & + + popcnt( iand( xor(key_in(i,1),iand(key_in(i,1),act_bitmask(i,1))), virt_bitmask(i,1) )) & + + popcnt( iand( xor(key_in(i,2),iand(key_in(i,2),act_bitmask(i,2))), virt_bitmask(i,2) )) enddo endif is_a_two_holes_two_particles = (i_diff >3) @@ -434,8 +220,8 @@ integer function number_of_holes_verbose(key_in) print*,'jey_in = ' call debug_det(key_in,N_int) number_of_holes_verbose = 0 - key_tmp(1,1) = xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))) - key_tmp(1,2) = xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,1,1))) + key_tmp(1,1) = xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))) + key_tmp(1,2) = xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,1))) call debug_det(key_tmp,N_int) key_tmp(1,1) = iand(key_tmp(1,1),reunion_of_core_inact_bitmask(1,1)) key_tmp(1,2) = iand(key_tmp(1,2),reunion_of_core_inact_bitmask(1,2)) @@ -446,8 +232,8 @@ integer function number_of_holes_verbose(key_in) ! number_of_holes_verbose = number_of_holes_verbose + popcnt(key_tmp(1,1)) & ! + popcnt(key_tmp(1,2)) number_of_holes_verbose = number_of_holes_verbose & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1)))), reunion_of_core_inact_bitmask(1,1)) )& + + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2)))), reunion_of_core_inact_bitmask(1,2)) ) print*,'----------------------' end @@ -464,8 +250,8 @@ integer function number_of_particles_verbose(key_in) print*,'jey_in = ' call debug_det(key_in,N_int) number_of_particles_verbose = 0 - key_tmp(1,1) = xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,1,1))) - key_tmp(1,2) = xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,1,1))) + key_tmp(1,1) = xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,1))) + key_tmp(1,2) = xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,1))) call debug_det(key_tmp,N_int) key_tmp(1,1) = iand(key_tmp(1,2),virt_bitmask(1,2)) key_tmp(1,2) = iand(key_tmp(1,2),virt_bitmask(1,2)) @@ -476,18 +262,16 @@ integer function number_of_particles_verbose(key_in) ! number_of_particles_verbose = number_of_particles_verbose + popcnt(key_tmp(1,1)) & ! + popcnt(key_tmp(1,2)) number_of_particles_verbose = number_of_particles_verbose & - + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & - + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) + + popcnt( iand( iand( xor(key_in(1,1),iand(key_in(1,1),act_bitmask(1,1))), virt_bitmask(1,1) ), virt_bitmask(1,1)) ) & + + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),act_bitmask(1,2))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) end logical function is_a_1h1p(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_1h1p = .False. - if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.1)then - is_a_1h1p = .True. - endif + + is_a_1h1p = (number_of_holes(key_in) == 1) .and. (number_of_particles(key_in) == 1) end @@ -495,10 +279,8 @@ logical function is_a_1h2p(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_1h2p = .False. - if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.2)then - is_a_1h2p = .True. - endif + + is_a_1h2p = (number_of_holes(key_in) == 1) .and. (number_of_particles(key_in) == 2) end @@ -506,10 +288,8 @@ logical function is_a_2h1p(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_2h1p = .False. - if(number_of_holes(key_in).eq.2 .and. number_of_particles(key_in).eq.1)then - is_a_2h1p = .True. - endif + + is_a_2h1p = (number_of_holes(key_in) == 2) .and. (number_of_particles(key_in) == 1) end @@ -517,10 +297,8 @@ logical function is_a_1h(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_1h = .False. - if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.0)then - is_a_1h = .True. - endif + + is_a_1h = (number_of_holes(key_in) == 1) .and. (number_of_particles(key_in) == 0) end @@ -528,10 +306,8 @@ logical function is_a_1p(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_1p = .False. - if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.1)then - is_a_1p = .True. - endif + + is_a_1p = (number_of_holes(key_in) == 0) .and. (number_of_particles(key_in) == 1) end @@ -539,10 +315,8 @@ logical function is_a_2p(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_2p = .False. - if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.2)then - is_a_2p = .True. - endif + + is_a_2p = (number_of_holes(key_in) == 0) .and. (number_of_particles(key_in) == 2) end @@ -550,10 +324,8 @@ logical function is_a_2h(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: number_of_particles, number_of_holes - is_a_2h = .False. - if(number_of_holes(key_in).eq.2 .and. number_of_particles(key_in).eq.0)then - is_a_2h = .True. - endif + + is_a_2h = (number_of_holes(key_in) == 2) .and. (number_of_particles(key_in) == 0) end diff --git a/src/bitmask/bitmasks.ezfio_config b/src/bitmask/bitmasks.ezfio_config index c133d8fe..dfb95c83 100644 --- a/src/bitmask/bitmasks.ezfio_config +++ b/src/bitmask/bitmasks.ezfio_config @@ -1,8 +1,4 @@ bitmasks N_int integer bit_kind integer - N_mask_gen integer - generators integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen) - N_mask_cas integer - cas integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,bitmasks_N_mask_cas) diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index bbcff63c..91617397 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -97,103 +97,9 @@ BEGIN_PROVIDER [ integer(bit_kind), ref_bitmask, (N_int,2)] ref_bitmask = HF_bitmask END_PROVIDER -BEGIN_PROVIDER [ integer, N_generators_bitmask ] - implicit none - BEGIN_DOC - ! Number of bitmasks for generators - END_DOC - logical :: exists - PROVIDE ezfio_filename N_int - - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_gen(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' - endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_generators_bitmask = 1 - endif - ASSERT (N_generators_bitmask > 0) - call write_int(6,N_generators_bitmask,'N_generators_bitmask') - 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( N_generators_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_generators_bitmask with MPI' - endif - IRP_ENDIF - - -END_PROVIDER -BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ] - implicit none - BEGIN_DOC - ! Number of bitmasks for generators - END_DOC - logical :: exists - PROVIDE ezfio_filename N_int - - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_gen(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' - endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_generators_bitmask_restart = 1 - endif - ASSERT (N_generators_bitmask_restart > 0) - call write_int(6,N_generators_bitmask_restart,'N_generators_bitmask_restart') - 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( N_generators_bitmask_restart, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_generators_bitmask_restart with MPI' - endif - IRP_ENDIF - - -END_PROVIDER - - - - -BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ] +BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6) ] implicit none BEGIN_DOC ! Bitmasks for generator determinants. @@ -215,231 +121,19 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_gen ! END_DOC logical :: exists - PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask N_int - PROVIDE generators_bitmask_restart + PROVIDE ezfio_filename full_ijkl_bitmask - if (mpi_master) then - call ezfio_has_bitmasks_generators(exists) - if (exists) then - call ezfio_get_bitmasks_generators(generators_bitmask_restart) - else - integer :: k, ispin - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask_restart(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,s_part ,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_part1,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_part2,k) = full_ijkl_bitmask(i) - enddo - enddo + integer :: ispin, i + do ispin=1,2 + do i=1,N_int + generators_bitmask(i,ispin,s_hole ) = reunion_of_inact_act_bitmask(i,ispin) + generators_bitmask(i,ispin,s_part ) = reunion_of_act_virt_bitmask(i,ispin) + generators_bitmask(i,ispin,d_hole1) = reunion_of_inact_act_bitmask(i,ispin) + generators_bitmask(i,ispin,d_part1) = reunion_of_act_virt_bitmask(i,ispin) + generators_bitmask(i,ispin,d_hole2) = reunion_of_inact_act_bitmask(i,ispin) + generators_bitmask(i,ispin,d_part2) = reunion_of_act_virt_bitmask(i,ispin) enddo - endif - - integer :: i - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask_restart(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_hole,k) ) - generators_bitmask_restart(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_part,k) ) - generators_bitmask_restart(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole1,k) ) - generators_bitmask_restart(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part1,k) ) - generators_bitmask_restart(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole2,k) ) - generators_bitmask_restart(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part2,k) ) - enddo - enddo - enddo - 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( generators_bitmask_restart, N_int*2*6*N_generators_bitmask_restart, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read generators_bitmask_restart with MPI' - endif - IRP_ENDIF - -END_PROVIDER - - -BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for generator determinants. - ! (N_int, alpha/beta, hole/particle, generator). - ! - ! 3rd index is : - ! - ! * 1 : hole for single exc - ! - ! * 2 : particle for single exc - ! - ! * 3 : hole for 1st exc of double - ! - ! * 4 : particle for 1st exc of double - ! - ! * 5 : hole for 2nd exc of double - ! - ! * 6 : particle for 2nd exc of double - ! - END_DOC - logical :: exists - PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask - - if (mpi_master) then - call ezfio_has_bitmasks_generators(exists) - if (exists) then - call ezfio_get_bitmasks_generators(generators_bitmask) - else - integer :: k, ispin, i - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,s_part ,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_part1,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_part2,k) = full_ijkl_bitmask(i) - enddo - enddo - enddo - endif - - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_hole,k) ) - generators_bitmask(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_part,k) ) - generators_bitmask(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole1,k) ) - generators_bitmask(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part1,k) ) - generators_bitmask(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole2,k) ) - generators_bitmask(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part2,k) ) - enddo - enddo - enddo - 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( generators_bitmask, N_int*2*6*N_generators_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read generators_bitmask with MPI' - endif - IRP_ENDIF - -END_PROVIDER - -BEGIN_PROVIDER [ integer, N_cas_bitmask ] - implicit none - BEGIN_DOC - ! Number of bitmasks for CAS - END_DOC - logical :: exists - PROVIDE ezfio_filename - PROVIDE N_cas_bitmask N_int - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_cas(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_cas(N_cas_bitmask) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' - endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_cas_bitmask = 1 - endif - call write_int(6,N_cas_bitmask,'N_cas_bitmask') - endif - ASSERT (N_cas_bitmask > 0) - 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( N_cas_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_cas_bitmask with MPI' - endif - IRP_ENDIF - -END_PROVIDER - -BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) - END_DOC - logical :: exists - integer :: i,i_part,i_gen,j,k - PROVIDE ezfio_filename generators_bitmask_restart full_ijkl_bitmask - PROVIDE n_generators_bitmask HF_bitmask - - if (mpi_master) then - call ezfio_has_bitmasks_cas(exists) - if (exists) then - call ezfio_get_bitmasks_cas(cas_bitmask) - else - if(N_generators_bitmask == 1)then - do j=1, N_cas_bitmask - do i=1, N_int - cas_bitmask(i,1,j) = iand(not(HF_bitmask(i,1)),full_ijkl_bitmask(i)) - cas_bitmask(i,2,j) = iand(not(HF_bitmask(i,2)),full_ijkl_bitmask(i)) - enddo - enddo - else - i_part = 2 - i_gen = 1 - do j=1, N_cas_bitmask - do i=1, N_int - cas_bitmask(i,1,j) = generators_bitmask_restart(i,1,i_part,i_gen) - cas_bitmask(i,2,j) = generators_bitmask_restart(i,2,i_part,i_gen) - enddo - enddo - endif - endif - do i=1,N_cas_bitmask - do j = 1, N_cas_bitmask - do k=1,N_int - cas_bitmask(k,j,i) = iand(cas_bitmask(k,j,i),full_ijkl_bitmask(k)) - enddo - enddo - enddo - write(*,*) 'Read CAS bitmask' - 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( cas_bitmask, N_int*2*N_cas_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read cas_bitmask with MPI' - endif - IRP_ENDIF - + enddo END_PROVIDER @@ -469,6 +163,19 @@ BEGIN_PROVIDER [integer(bit_kind), reunion_of_inact_act_bitmask, (N_int,2)] enddo END_PROVIDER +BEGIN_PROVIDER [integer(bit_kind), reunion_of_act_virt_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive and active bitmasks + END_DOC + integer :: i,j + + do i = 1, N_int + reunion_of_act_virt_bitmask(i,1) = ior(virt_bitmask(i,1),act_bitmask(i,1)) + reunion_of_act_virt_bitmask(i,2) = ior(virt_bitmask(i,2),act_bitmask(i,2)) + enddo +END_PROVIDER + BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] implicit none @@ -491,8 +198,8 @@ BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] END_DOC integer :: i,j do i = 1, N_int - reunion_of_bitmask(i,1) = ior(ior(cas_bitmask(i,1,1),inact_bitmask(i,1)),virt_bitmask(i,1)) - reunion_of_bitmask(i,2) = ior(ior(cas_bitmask(i,2,1),inact_bitmask(i,2)),virt_bitmask(i,2)) + reunion_of_bitmask(i,1) = ior(ior(act_bitmask(i,1),inact_bitmask(i,1)),virt_bitmask(i,1)) + reunion_of_bitmask(i,2) = ior(ior(act_bitmask(i,2),inact_bitmask(i,2)),virt_bitmask(i,2)) enddo END_PROVIDER @@ -512,14 +219,6 @@ END_PROVIDER enddo END_PROVIDER -BEGIN_PROVIDER [ integer, i_bitmask_gen ] - implicit none - BEGIN_DOC - ! Current bitmask for the generators - END_DOC - i_bitmask_gen = 1 -END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), unpaired_alpha_electrons, (N_int)] implicit none @@ -537,21 +236,7 @@ BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] implicit none integer :: i,j do i = 1, N_int - closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),cas_bitmask(i,1,1)) - closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),cas_bitmask(i,2,1)) + closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),act_bitmask(i,1)) + closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),act_bitmask(i,2)) enddo END_PROVIDER - - -BEGIN_PROVIDER [ integer(bit_kind), reunion_of_cas_inact_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks - END_DOC - integer :: i,j - do i = 1, N_int - reunion_of_cas_inact_bitmask(i,1) = ior(act_bitmask(i,1),inact_bitmask(i,1)) - reunion_of_cas_inact_bitmask(i,2) = ior(act_bitmask(i,2),inact_bitmask(i,2)) - enddo -END_PROVIDER - diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index b016f1fd..d30e989f 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -129,6 +129,15 @@ BEGIN_PROVIDER [integer, dim_list_inact_orb] dim_list_inact_orb = max(n_inact_orb,1) END_PROVIDER +BEGIN_PROVIDER [integer, dim_list_core_inact_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_core. + ! it is at least 1 + END_DOC + dim_list_core_inact_orb = max(n_core_inact_orb,1) +END_PROVIDER + BEGIN_PROVIDER [integer, dim_list_act_orb] implicit none BEGIN_DOC @@ -168,43 +177,67 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), core_bitmask , (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), act_bitmask , (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask , (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), del_bitmask , (N_int,2) ] implicit none BEGIN_DOC - ! Bitmask identifying the core/inactive/active/virtual/deleted MOs + ! Bitmask identifying the core MOs END_DOC - core_bitmask = 0_bit_kind - inact_bitmask = 0_bit_kind - act_bitmask = 0_bit_kind - virt_bitmask = 0_bit_kind - del_bitmask = 0_bit_kind - if(n_core_orb > 0)then call list_to_bitstring( core_bitmask(1,1), list_core, n_core_orb, N_int) call list_to_bitstring( core_bitmask(1,2), list_core, n_core_orb, N_int) endif + END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the inactive MOs + END_DOC + inact_bitmask = 0_bit_kind if(n_inact_orb > 0)then call list_to_bitstring( inact_bitmask(1,1), list_inact, n_inact_orb, N_int) call list_to_bitstring( inact_bitmask(1,2), list_inact, n_inact_orb, N_int) endif + END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), act_bitmask , (N_int,2) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the active MOs + END_DOC + act_bitmask = 0_bit_kind if(n_act_orb > 0)then call list_to_bitstring( act_bitmask(1,1), list_act, n_act_orb, N_int) call list_to_bitstring( act_bitmask(1,2), list_act, n_act_orb, N_int) endif + END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask , (N_int,2) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the virtual MOs + END_DOC + virt_bitmask = 0_bit_kind if(n_virt_orb > 0)then call list_to_bitstring( virt_bitmask(1,1), list_virt, n_virt_orb, N_int) call list_to_bitstring( virt_bitmask(1,2), list_virt, n_virt_orb, N_int) endif + END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), del_bitmask , (N_int,2) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the deleted MOs + END_DOC + + del_bitmask = 0_bit_kind + if(n_del_orb > 0)then call list_to_bitstring( del_bitmask(1,1), list_del, n_del_orb, N_int) call list_to_bitstring( del_bitmask(1,2), list_del, n_del_orb, N_int) endif -END_PROVIDER + END_PROVIDER @@ -322,13 +355,12 @@ END_PROVIDER enddo print *, 'Active MOs:' print *, list_act(1:n_act_orb) - print*, list_act_reverse(1:n_act_orb) END_PROVIDER - BEGIN_PROVIDER [ integer, list_core_inact , (n_core_inact_orb) ] + BEGIN_PROVIDER [ integer, list_core_inact , (dim_list_core_inact_orb) ] &BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_num) ] implicit none BEGIN_DOC diff --git a/src/bitmask/modify_bitmasks.irp.f b/src/bitmask/modify_bitmasks.irp.f index fa660680..834be6c8 100644 --- a/src/bitmask/modify_bitmasks.irp.f +++ b/src/bitmask/modify_bitmasks.irp.f @@ -1,26 +1,5 @@ use bitmasks -subroutine initialize_bitmask_to_restart_ones - implicit none - integer :: i,j,k,l,m - integer :: ispin - BEGIN_DOC - ! Initialization of the generators_bitmask to the restart bitmask - END_DOC - do i = 1, N_int - do k=1,N_generators_bitmask - do ispin=1,2 - generators_bitmask(i,ispin,s_hole ,k) = generators_bitmask_restart(i,ispin,s_hole ,k) - generators_bitmask(i,ispin,s_part ,k) = generators_bitmask_restart(i,ispin,s_part ,k) - generators_bitmask(i,ispin,d_hole1,k) = generators_bitmask_restart(i,ispin,d_hole1,k) - generators_bitmask(i,ispin,d_part1,k) = generators_bitmask_restart(i,ispin,d_part1,k) - generators_bitmask(i,ispin,d_hole2,k) = generators_bitmask_restart(i,ispin,d_hole2,k) - generators_bitmask(i,ispin,d_part2,k) = generators_bitmask_restart(i,ispin,d_part2,k) - enddo - enddo - enddo -end - subroutine modify_bitmasks_for_hole(i_hole) implicit none @@ -33,26 +12,22 @@ subroutine modify_bitmasks_for_hole(i_hole) END_DOC ! Set to Zero the holes - do k=1,N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_holes_bitmask(l) do ispin=1,2 do j = 1, N_int - generators_bitmask(j,ispin,i,k) = 0_bit_kind + generators_bitmask(j,ispin,i) = 0_bit_kind enddo enddo - enddo enddo k = shiftr(i_hole-1,bit_kind_shift)+1 j = i_hole-shiftl(k-1,bit_kind_shift)-1 - do m = 1, N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_holes_bitmask(l) do ispin=1,2 - generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + generators_bitmask(k,ispin,i) = ibset(generators_bitmask(k,ispin,i),j) enddo - enddo enddo end @@ -69,13 +44,11 @@ subroutine modify_bitmasks_for_hole_in_out(i_hole) k = shiftr(i_hole-1,bit_kind_shift)+1 j = i_hole-shiftl(k-1,bit_kind_shift)-1 - do m = 1, N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_holes_bitmask(l) do ispin=1,2 - generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + generators_bitmask(k,ispin,i) = ibset(generators_bitmask(k,ispin,i),j) enddo - enddo enddo end @@ -91,75 +64,67 @@ subroutine modify_bitmasks_for_particl(i_part) END_DOC ! Set to Zero the particles - do k=1,N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_particl_bitmask(l) - do ispin=1,2 + do ispin=1,2 do j = 1, N_int - generators_bitmask(j,ispin,i,k) = 0_bit_kind + generators_bitmask(j,ispin,i) = 0_bit_kind enddo - enddo enddo enddo k = shiftr(i_part-1,bit_kind_shift)+1 j = i_part-shiftl(k-1,bit_kind_shift)-1 - do m = 1, N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_particl_bitmask(l) do ispin=1,2 - generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + generators_bitmask(k,ispin,i) = ibset(generators_bitmask(k,ispin,i),j) enddo - enddo enddo end -subroutine set_bitmask_particl_as_input(input_bimask) +subroutine set_bitmask_particl_as_input(input_bitmask) implicit none - integer(bit_kind), intent(in) :: input_bimask(N_int,2) + integer(bit_kind), intent(in) :: input_bitmask(N_int,2) integer :: i,j,k,l,m integer :: ispin BEGIN_DOC ! set the generators_bitmask for the particles -! as the input_bimask +! as the input_bitmask END_DOC - do k=1,N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_particl_bitmask(l) - do ispin=1,2 + do ispin=1,2 do j = 1, N_int - generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin) + generators_bitmask(j,ispin,i) = input_bitmask(j,ispin) enddo enddo - enddo enddo touch generators_bitmask end -subroutine set_bitmask_hole_as_input(input_bimask) +subroutine set_bitmask_hole_as_input(input_bitmask) implicit none - integer(bit_kind), intent(in) :: input_bimask(N_int,2) + integer(bit_kind), intent(in) :: input_bitmask(N_int,2) integer :: i,j,k,l,m integer :: ispin BEGIN_DOC ! set the generators_bitmask for the holes -! as the input_bimask +! as the input_bitmask END_DOC - do k=1,N_generators_bitmask - do l = 1, 3 + do l = 1, 3 i = index_holes_bitmask(l) do ispin=1,2 do j = 1, N_int - generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin) + generators_bitmask(j,ispin,i) = input_bitmask(j,ispin) enddo enddo - enddo enddo touch generators_bitmask @@ -173,11 +138,10 @@ subroutine print_generators_bitmasks_holes allocate(key_tmp(N_int,2)) do l = 1, 3 - k = 1 - i = index_holes_bitmask(l) + i = index_holes_bitmask(l) do j = 1, N_int - key_tmp(j,1) = generators_bitmask(j,1,i,k) - key_tmp(j,2) = generators_bitmask(j,2,i,k) + key_tmp(j,1) = generators_bitmask(j,1,i) + key_tmp(j,2) = generators_bitmask(j,2,i) enddo print*,'' print*,'index hole = ',i @@ -195,57 +159,10 @@ subroutine print_generators_bitmasks_particles allocate(key_tmp(N_int,2)) do l = 1, 3 - k = 1 - i = index_particl_bitmask(l) + i = index_particl_bitmask(l) do j = 1, N_int - key_tmp(j,1) = generators_bitmask(j,1,i,k) - key_tmp(j,2) = generators_bitmask(j,2,i,k) - enddo - print*,'' - print*,'index particl ',i - call print_det(key_tmp,N_int) - print*,'' - enddo - deallocate(key_tmp) - -end - -subroutine print_generators_bitmasks_holes_for_one_generator(i_gen) - implicit none - integer, intent(in) :: i_gen - integer :: i,j,k,l - integer(bit_kind),allocatable :: key_tmp(:,:) - - allocate(key_tmp(N_int,2)) - do l = 1, 3 - k = i_gen - i = index_holes_bitmask(l) - do j = 1, N_int - key_tmp(j,1) = generators_bitmask(j,1,i,k) - key_tmp(j,2) = generators_bitmask(j,2,i,k) - enddo - print*,'' - print*,'index hole = ',i - call print_det(key_tmp,N_int) - print*,'' - enddo - deallocate(key_tmp) - -end - -subroutine print_generators_bitmasks_particles_for_one_generator(i_gen) - implicit none - integer, intent(in) :: i_gen - integer :: i,j,k,l - integer(bit_kind),allocatable :: key_tmp(:,:) - - allocate(key_tmp(N_int,2)) - do l = 1, 3 - k = i_gen - i = index_particl_bitmask(l) - do j = 1, N_int - key_tmp(j,1) = generators_bitmask(j,1,i,k) - key_tmp(j,2) = generators_bitmask(j,2,i,k) + key_tmp(j,1) = generators_bitmask(j,1,i) + key_tmp(j,2) = generators_bitmask(j,2,i) enddo print*,'' print*,'index particl ',i @@ -257,7 +174,7 @@ subroutine print_generators_bitmasks_particles_for_one_generator(i_gen) end - BEGIN_PROVIDER [integer, index_holes_bitmask, (3)] +BEGIN_PROVIDER [integer, index_holes_bitmask, (3)] implicit none BEGIN_DOC ! Index of the holes in the generators_bitmasks diff --git a/src/casscf/50.casscf.bats b/src/casscf/50.casscf.bats new file mode 100644 index 00000000..a0db725d --- /dev/null +++ b/src/casscf/50.casscf.bats @@ -0,0 +1,49 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh +source $QP_ROOT/quantum_package.rc + + +function run_stoch() { + thresh=$2 + test_exe casscf || skip + qp set perturbation do_pt2 True + qp set determinants n_det_max $3 + qp set davidson threshold_davidson 1.e-10 + qp set davidson n_states_diag 4 + qp run casscf | tee casscf.out + energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)" + eq $energy1 $1 $thresh +} + +@test "F2" { # 18.0198s + rm -rf f2_casscf + qp_create_ezfio -b aug-cc-pvdz ../input/f2.zmt -o f2_casscf + qp set_file f2_casscf + qp run scf + qp set_mo_class --core="[1-6,8-9]" --act="[7,10]" --virt="[11-46]" + run_stoch -198.773366970 1.e-4 100000 +} + +@test "N2" { # 18.0198s + rm -rf n2_casscf + qp_create_ezfio -b aug-cc-pvdz ../input/n2.xyz -o n2_casscf + qp set_file n2_casscf + qp run scf + qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" + run_stoch -109.0961643162 1.e-4 100000 +} + +@test "N2_stretched" { + rm -rf n2_stretched_casscf + qp_create_ezfio -b aug-cc-pvdz -m 7 ../input/n2_stretched.xyz -o n2_stretched_casscf + qp set_file n2_stretched_casscf + qp run scf | tee scf.out + qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" + qp set electrons elec_alpha_num 7 + qp set electrons elec_beta_num 7 + run_stoch -108.7860471300 1.e-4 100000 +# + +} + diff --git a/src/casscf/EZFIO.cfg b/src/casscf/EZFIO.cfg index ce51a064..4e4d3d3a 100644 --- a/src/casscf/EZFIO.cfg +++ b/src/casscf/EZFIO.cfg @@ -16,4 +16,16 @@ doc: If true, the CASSCF starts with a CISD wave function interface: ezfio,provider,ocaml default: True +[state_following_casscf] +type: logical +doc: If |true|, the CASSCF will try to follow the guess CI vector and orbitals +interface: ezfio,provider,ocaml +default: False + + +[level_shift_casscf] +type: Positive_float +doc: Energy shift on the virtual MOs to improve SCF convergence +interface: ezfio,provider,ocaml +default: 0.005 diff --git a/src/casscf/MORALITY b/src/casscf/MORALITY new file mode 100644 index 00000000..9701a647 --- /dev/null +++ b/src/casscf/MORALITY @@ -0,0 +1 @@ +the CASCF can be obtained if a proper guess is given to the WF part diff --git a/src/casscf/NEED b/src/casscf/NEED index b992ff71..d9da718e 100644 --- a/src/casscf/NEED +++ b/src/casscf/NEED @@ -1,4 +1,4 @@ cipsi selectors_full -generators_fluid +generators_cas two_body_rdm diff --git a/src/casscf/bavard.irp.f b/src/casscf/bavard.irp.f index 402e67ec..463c3ea4 100644 --- a/src/casscf/bavard.irp.f +++ b/src/casscf/bavard.irp.f @@ -1,6 +1,6 @@ ! -*- F90 -*- BEGIN_PROVIDER [logical, bavard] -! bavard=.true. - bavard=.false. +! bavard=.true. + bavard=.false. END_PROVIDER diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 8b5a365f..d83aa271 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -3,99 +3,27 @@ program casscf BEGIN_DOC ! TODO : Put the documentation of the program here END_DOC + call reorder_orbitals_for_casscf no_vvvv_integrals = .True. - SOFT_TOUCH no_vvvv_integrals - threshold_davidson = 1.d-7 - touch threshold_davidson - if(cisd_guess)then - logical :: converged - integer :: iteration - double precision :: energy - print*,'*******************************' - print*,'*******************************' - print*,'*******************************' - print*,'USING A CISD WAVE FUNCTION AS GUESS FOR THE MCSCF WF' - print*,'*******************************' - print*,'*******************************' - converged = .False. - iteration = 0 - generators_type = "HF" - touch generators_type - read_wf = .False. - touch read_wf - logical :: do_cisdtq - do_cisdtq = .True. - double precision :: thr - thr = 5.d-3 - do while (.not.converged) - call cisd_scf_iteration(converged,iteration,energy,thr) - if(HF_index.ne.1.and.iteration.gt.0)then - print*,'*******************************' - print*,'*******************************' - print*,'The HF determinant is not the dominant determinant in the CISD WF ...' - print*,'Therefore we skip the CISD WF ..' - print*,'*******************************' - print*,'*******************************' - do_cisdtq = .False. - exit - endif - if(iteration.gt.15.and..not.converged)then - print*,'It seems that the orbital optimization for the CISD WAVE FUNCTION CANNOT CONVERGE ...' - print*,'Passing to CISDTQ WAVE FUNCTION' - exit - endif - enddo - if(do_cisdtq)then - print*,'*******************************' - print*,'*******************************' - print*,'*******************************' - print*,'SWITCHING WITH A CISDTQ WAVE FUNCTION AS GUESS FOR THE MCSCF WF' - print*,'*******************************' - print*,'*******************************' - converged = .False. - iteration = 0 - read_wf = .False. - touch read_wf - pt2_max = 0.01d0 - touch pt2_max - energy = 0.d0 - do while (.not.converged) - call cisdtq_scf_iteration(converged,iteration,energy,thr) - if(HF_index.ne.1.and.iteration.gt.0)then - print*,'*******************************' - print*,'*******************************' - print*,'The HF determinant is not the dominant determinant in the CISDTQ WF ...' - print*,'Therefore we skip the CISDTQ WF ..' - print*,'*******************************' - print*,'*******************************' - exit - endif - if(iteration.gt.15.and..not.converged)then - print*,'It seems that the orbital optimization for the CISDTQ WAVE FUNCTION CANNOT CONVERGE ...' - print*,'Passing to CISDTQ WAVE FUNCTION' - exit - endif - enddo - endif - endif - read_wf = .False. - touch read_wf - pt2_max = 0.0d0 - touch pt2_max -! call run_cipsi_scf - call run + pt2_max = 0.02 + SOFT_TOUCH no_vvvv_integrals pt2_max + call run_stochastic_cipsi + call run end subroutine run implicit none double precision :: energy_old, energy - logical :: converged + logical :: converged,state_following_casscf_save integer :: iteration converged = .False. energy = 0.d0 mo_label = "MCSCF" iteration = 1 + state_following_casscf_save = state_following_casscf + state_following_casscf = .True. + touch state_following_casscf do while (.not.converged) call run_stochastic_cipsi energy_old = energy @@ -107,17 +35,22 @@ subroutine run call write_double(6,energy_improvement, 'Predicted energy improvement') converged = dabs(energy_improvement) < thresh_scf -! pt2_max = dabs(energy_improvement / pt2_relative_error) + pt2_max = dabs(energy_improvement / pt2_relative_error) mo_coef = NewOrbs + mo_occ = occnum call save_mos iteration += 1 - N_det = N_det/2 + N_det = max(N_det/2 ,N_states) psi_det = psi_det_sorted psi_coef = psi_coef_sorted read_wf = .True. call clear_mo_map SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef + if(iteration .gt. 3)then + state_following_casscf = state_following_casscf_save + touch state_following_casscf + endif enddo diff --git a/src/casscf/change_bitmasks.irp.f b/src/casscf/change_bitmasks.irp.f deleted file mode 100644 index cad6ec38..00000000 --- a/src/casscf/change_bitmasks.irp.f +++ /dev/null @@ -1,14 +0,0 @@ -subroutine only_act_bitmask - implicit none - integer :: i,j,k - do k = 1, N_generators_bitmask - do j = 1, 6 - do i = 1, N_int - generators_bitmask(i,1,j,k) = act_bitmask(i,1) - generators_bitmask(i,2,j,k) = act_bitmask(i,2) - enddo - enddo - enddo - touch generators_bitmask -end - diff --git a/src/casscf/cipsi_routines.irp.f b/src/casscf/cipsi_routines.irp.f deleted file mode 100644 index 272a7116..00000000 --- a/src/casscf/cipsi_routines.irp.f +++ /dev/null @@ -1,75 +0,0 @@ -subroutine run_cipsi_scf - implicit none - double precision :: energy_old, energy, extrap,extrap_old,pt2_max_begin - logical :: converged - integer :: iteration - print*,'*********************************' - print*,'*********************************' - print*,' DOING THE CIPSI-SCF ' - print*,'*********************************' - converged = .False. - pt2_max_begin = pt2_max - energy = 0.d0 - extrap = 0.d0 - mo_label = "MCSCF" - iteration = 1 - threshold_davidson = 1.d-09 - touch threshold_davidson - do while (.not.converged) - print*,'' - call write_int(6,iteration,'CI STEP OF THE ITERATION = ') - call write_double(6,pt2_max,'PT2 MAX = ') - !call cisd_guess_wf - generators_type = "CAS" - touch generators_type - call run_stochastic_cipsi - call change_orb_cipsi(converged,iteration,energy) - if(iteration.gt.n_it_scf_max.and..not.converged)then - print*,'It seems that the orbital optimization for the CISDTQ WAVE FUNCTION CANNOT CONVERGE ...' - print*,'The required delta E was :',thresh_scf - print*,'The obtained delta E was :',extrap - extrap_old - print*,'After ',iteration,'iterations ...' - print*,'Getting out of the SCF loop ...' - exit - endif - iteration += 1 - enddo - -end - -subroutine change_orb_cipsi(converged,iteration,energy) - implicit none - double precision :: energy_old, extrap,extrap_old,pt2_max_begin - double precision, intent(inout):: energy - logical, intent(out) :: converged - integer, intent(in) :: iteration - extrap_old = energy - energy = eone+etwo+ecore - extrap = extrapolated_energy(2,1) - - call write_time(6) - call write_int(6,iteration,'CAS-SCF iteration') - call write_double(6,energy,'CAS-SCF variational energy') - call write_double(6,extrap,'CAS-SCF extrapolated energy') - call write_double(6,extrap - extrap_old,'Change in extrapolated energy') - energy = extrap - call write_double(6,energy_improvement, 'Predicted energy improvement') - - converged = dabs(extrap - extrap_old) < thresh_scf - pt2_max = dabs(extrap - extrap_old) * 10.d0 - pt2_max = min(pt2_max,1.d-2) - pt2_max = max(pt2_max,1.d-10) - if(N_det.gt.10**6)then - pt2_max = max(pt2_max,1.d-2) - endif - - mo_coef = NewOrbs - call save_mos - call map_deinit(mo_integrals_map) - N_det = N_det/2 - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - read_wf = .True. - FREE mo_integrals_map mo_two_e_integrals_in_map - SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef -end diff --git a/src/casscf/cisd_routine.irp.f b/src/casscf/cisd_routine.irp.f deleted file mode 100644 index a4cbfcfb..00000000 --- a/src/casscf/cisd_routine.irp.f +++ /dev/null @@ -1,85 +0,0 @@ -subroutine cisd_scf_iteration(converged,iteration,energy,thr) - implicit none - double precision, intent(in) :: thr - logical, intent(out) :: converged - integer, intent(inout) :: iteration - double precision, intent(out) :: energy - converged = .False. - call only_act_bitmask - N_det = N_det_generators - psi_coef = psi_coef_generators - psi_det = psi_det_generators - touch N_det psi_coef psi_det - call run_cisd - call change_orb_cisd(converged,iteration,energy,thr) -end - - -subroutine cisd_guess_wf - implicit none - call only_act_bitmask - N_det = N_det_generators - psi_coef = psi_coef_generators - psi_det = psi_det_generators - touch N_det psi_coef psi_det - generators_type = "HF" - touch generators_type - call run_cisd - touch N_det psi_coef psi_det psi_coef_sorted psi_det_sorted psi_det_sorted_order psi_average_norm_contrib_sorted - -end - - - -subroutine change_orb_cisd(converged,iteration,energy,thr) - implicit none - double precision, intent(in) :: thr - logical, intent(inout) :: converged - integer, intent(inout) :: iteration - double precision, intent(inout) :: energy - double precision :: energy_old - energy_old = energy - - energy = eone+etwo+ecore - - call write_time(6) - call write_int(6,iteration,'CISD-SCF iteration') - call write_double(6,energy,'CISD-SCF energy') - call write_double(6,energy_improvement, 'Predicted energy improvement') - converged = dabs(energy_improvement) < thr - - mo_coef = NewOrbs - call save_mos - call map_deinit(mo_integrals_map) - FREE mo_integrals_map mo_two_e_integrals_in_map - iteration += 1 - -end - -subroutine run_cisd - implicit none - integer :: i - - if(pseudo_sym)then - call H_apply_cisd_sym - else - call H_apply_cisd - endif - print *, 'N_det = ', N_det - print*,'******************************' - print *, 'Energies of the states:' - do i = 1,N_states - print *, i, CI_energy(i) - enddo - if (N_states > 1) then - print*,'******************************' - print*,'Excitation energies ' - do i = 2, N_states - print*, i ,CI_energy(i) - CI_energy(1) - enddo - endif - psi_coef = ci_eigenvectors - SOFT_TOUCH psi_coef - call save_wavefunction - -end diff --git a/src/casscf/cisdtq_routine.irp.f b/src/casscf/cisdtq_routine.irp.f deleted file mode 100644 index 0479d462..00000000 --- a/src/casscf/cisdtq_routine.irp.f +++ /dev/null @@ -1,47 +0,0 @@ -subroutine cisdtq_scf_iteration(converged,iteration,energy,thr) - implicit none - double precision, intent(in) :: thr - logical, intent(out) :: converged - integer, intent(inout) :: iteration - double precision, intent(inout) :: energy - converged = .False. - call only_act_bitmask - generators_type = "HF_SD" - threshold_generators = 0.99d0 - touch threshold_generators - touch generators_type - selection_factor = 5 - touch selection_factor - call run_stochastic_cipsi - call change_orb_cisdtq(converged,iteration,energy,thr) -end - -subroutine change_orb_cisdtq(converged,iteration,energy,thr) - implicit none - double precision, intent(in) :: thr - logical, intent(inout) :: converged - integer, intent(inout) :: iteration - double precision, intent(inout) :: energy - double precision :: extrap,extrap_old,pt2_max_begin - extrap_old = energy - extrap = extrapolated_energy(2,1) - energy = extrap - - call write_time(6) - call write_int(6,iteration,'CISDTQ-SCF iteration') - call write_double(6,energy,'CISDTQ-SCF variational energy') - call write_double(6,extrap,'CISDTQ-SCF extrapolated energy') - call write_double(6,extrap - extrap_old,'Change in extrapolated energy') - - converged = dabs(extrap - extrap_old) < thr - pt2_max = dabs(extrap - extrap_old) * 10.d0 - pt2_max = max(pt2_max,1.d-10) - - mo_coef = NewOrbs - call save_mos - call map_deinit(mo_integrals_map) - FREE mo_integrals_map mo_two_e_integrals_in_map - iteration += 1 - -end - diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 292067b4..3d1ff0f9 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -56,8 +56,8 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] uu = list_act(u) do t = 1, n_act_orb tt = list_act(t) -! P0tuvx(t,u,v,x) = state_av_act_two_rdm_openmp_spin_trace_mo(t,v,u,x) P0tuvx(t,u,v,x) = state_av_act_two_rdm_spin_trace_mo(t,v,u,x) +! P0tuvx(t,u,v,x) = act_two_rdm_spin_trace_mo(t,v,u,x) enddo enddo enddo diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 8dce87aa..362da85d 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -5,10 +5,57 @@ program print_2rdm ! ! useful to test the active part of the spin trace 2 rdms END_DOC - no_vvvv_integrals = .True. +!no_vvvv_integrals = .True. read_wf = .True. - touch read_wf no_vvvv_integrals - call routine +!touch read_wf no_vvvv_integrals +!call routine +!call routine_bis + call print_grad +end + +subroutine print_grad + implicit none + integer :: i + do i = 1, nMonoEx + if(dabs(gradvec2(i)).gt.1.d-5)then + print*,'' + print*,i,gradvec2(i),excit(:,i) + endif + enddo +end + +subroutine routine_bis + implicit none + integer :: i,j + double precision :: accu_d,accu_od +!accu_d = 0.d0 +!accu_od = 0.d0 +!print*,'' +!print*,'' +!print*,'' +!do i = 1, mo_num +! write(*,'(100(F8.5,X))')super_ci_dm(i,:) +! accu_d += super_ci_dm(i,i) +! do j = i+1, mo_num +! accu_od += dabs(super_ci_dm(i,j) - super_ci_dm(j,i)) +! enddo +!enddo +!print*,'' +!print*,'' +!print*,'accu_d = ',accu_d +!print*,'n_elec = ',elec_num +!print*,'accu_od= ',accu_od +!print*,'' +!accu_d = 0.d0 +!do i = 1, N_det +! accu_d += psi_coef(i,1)**2 +!enddo +!print*,'accu_d = ',accu_d +!provide superci_natorb + + provide switch_mo_coef + mo_coef = switch_mo_coef + call save_mos end subroutine routine diff --git a/src/casscf/grad_old.irp.f b/src/casscf/grad_old.irp.f new file mode 100644 index 00000000..d60a60c8 --- /dev/null +++ b/src/casscf/grad_old.irp.f @@ -0,0 +1,74 @@ + +BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)] + BEGIN_DOC + ! calculate the orbital gradient by hand, i.e. for + ! each determinant I we determine the string E_pq |I> (alpha and beta + ! separately) and generate + ! sum_I c_I is then the pq component of the orbital + ! gradient + ! E_pq = a^+_pa_q + a^+_Pa_Q + END_DOC + implicit none + integer :: ii,tt,aa,indx,ihole,ipart,istate + real*8 :: res + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + call calc_grad_elem(ihole,ipart,res) + gradvec_old(indx)=res + end do + + real*8 :: norm_grad + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec_old(indx)*gradvec_old(indx) + end do + norm_grad=sqrt(norm_grad) + if (bavard) then + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + endif + + +END_PROVIDER + +subroutine calc_grad_elem(ihole,ipart,res) + BEGIN_DOC + ! eq 18 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 , q=hole, p=particle + END_DOC + implicit none + integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) + real*8 :: i_H_psi_array(N_states),phase + allocate(det_mu(N_int,2)) + allocate(det_mu_ex(N_int,2)) + + res=0.D0 + + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation on it + call det_copy(det_mu,det_mu_ex,N_int) + call do_signed_mono_excitation(det_mu,det_mu_ex,nu & + ,ihole,ipart,ispin,phase,ierr) + if (ierr.eq.1) then + call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase + end do + end if + end do + end do + + ! state-averaged gradient + res*=2.D0/dble(N_states) + +end subroutine calc_grad_elem + diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f index f00bc7c8..e717e822 100644 --- a/src/casscf/gradient.irp.f +++ b/src/casscf/gradient.irp.f @@ -60,79 +60,6 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)] - BEGIN_DOC - ! calculate the orbital gradient by hand, i.e. for - ! each determinant I we determine the string E_pq |I> (alpha and beta - ! separately) and generate - ! sum_I c_I is then the pq component of the orbital - ! gradient - ! E_pq = a^+_pa_q + a^+_Pa_Q - END_DOC - implicit none - integer :: ii,tt,aa,indx,ihole,ipart,istate - real*8 :: res - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - call calc_grad_elem(ihole,ipart,res) - gradvec(indx)=res - end do - - real*8 :: norm_grad - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec(indx)*gradvec(indx) - end do - norm_grad=sqrt(norm_grad) - if (bavard) then - write(6,*) - write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad - write(6,*) - endif - - -END_PROVIDER - -subroutine calc_grad_elem(ihole,ipart,res) - BEGIN_DOC - ! eq 18 of Siegbahn et al, Physica Scripta 1980 - ! we calculate 2 , q=hole, p=particle - END_DOC - implicit none - integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) - real*8 :: i_H_psi_array(N_states),phase - allocate(det_mu(N_int,2)) - allocate(det_mu_ex(N_int,2)) - - res=0.D0 - - do mu=1,n_det - ! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 - ! do the monoexcitation on it - call det_copy(det_mu,det_mu_ex,N_int) - call do_signed_mono_excitation(det_mu,det_mu_ex,nu & - ,ihole,ipart,ispin,phase,ierr) - if (ierr.eq.1) then - call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase - end do - end if - end do - end do - - ! state-averaged gradient - res*=2.D0/dble(N_states) - -end subroutine calc_grad_elem - BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] BEGIN_DOC ! calculate the orbital gradient from density @@ -171,11 +98,9 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] norm_grad+=gradvec2(indx)*gradvec2(indx) end do norm_grad=sqrt(norm_grad) -! if (bavard) then write(6,*) write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad write(6,*) -! endif END_PROVIDER diff --git a/src/casscf/h_apply.irp.f b/src/casscf/h_apply.irp.f deleted file mode 100644 index 6fcb2900..00000000 --- a/src/casscf/h_apply.irp.f +++ /dev/null @@ -1,18 +0,0 @@ -! Generates subroutine H_apply_cisd -! ---------------------------------- - -BEGIN_SHELL [ /usr/bin/env python2 ] -from generate_h_apply import H_apply -H = H_apply("cisd",do_double_exc=True) -print H - -from generate_h_apply import H_apply -H = H_apply("cisdtq",do_double_exc=True) -H.set_selection_pt2("epstein_nesbet_2x2") -print H - -H = H_apply("cisd_sym",do_double_exc=True) -H.filter_only_connected_to_hf() -print H -END_SHELL - diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f index 06aed6ef..52be1b76 100644 --- a/src/casscf/hessian.irp.f +++ b/src/casscf/hessian.irp.f @@ -536,9 +536,6 @@ real*8 function hessmat_taub(t,a,u,b) integer :: v3,x3 real*8 :: term,t1,t2,t3 - double precision,allocatable :: P0tuvx_no_t(:,:,:) - double precision :: bielec_pqxx_no_2(n_act_orb,n_act_orb) - double precision :: bielec_pxxq_no_2(n_act_orb,n_act_orb) tt=list_act(t) aa=list_virt(a) if (t == u) then @@ -548,87 +545,59 @@ real*8 function hessmat_taub(t,a,u,b) t2=0.D0 t3=0.D0 t1-=occnum(tt)*Fipq(tt,tt) - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - t2+=P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) - end do - end do do v=1,n_act_orb vv=list_act(v) v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) x3=x+n_core_inact_orb - t2+=(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq_no(aa,x3,v3,aa) - end do - end do - do y=1,n_act_orb - do x=1,n_act_orb - xx=list_act(x) - do v=1,n_act_orb - t3-=P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) + t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq_no(aa,x3,v3,aa)) + do y=1,n_act_orb + t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) end do end do end do - term=t1+2.d0*(t2+t3) + term=t1+t2+t3 else bb=list_virt(b) ! ta/tb b/=a - term=0.5d0*occnum(tt)*Fipq(aa,bb) - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - term = term + P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) - end do - end do + term=occnum(tt)*Fipq(aa,bb) do v=1,n_act_orb vv=list_act(v) v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) x3=x+n_core_inact_orb - term= term + (P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & - *bielec_pxxq_no(aa,x3,v3,bb) + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & + *bielec_pxxq_no(aa,x3,v3,bb)) end do end do - term += term end if else ! ta/ub t/=u uu=list_act(u) bb=list_virt(b) - allocate(P0tuvx_no_t(n_act_orb,n_act_orb,n_act_orb)) - P0tuvx_no_t(:,:,:) = P0tuvx_no(t,:,:,:) - do x=1,n_act_orb - x3=x+n_core_inact_orb - do v=1,n_act_orb - v3=v+n_core_inact_orb - bielec_pqxx_no_2(v,x) = bielec_pqxx_no(aa,bb,v3,x3) - bielec_pxxq_no_2(v,x) = bielec_pxxq_no(aa,v3,x3,bb) - end do - end do term=0.D0 - do x=1,n_act_orb - do v=1,n_act_orb - term += P0tuvx_no_t(u,v,x)*bielec_pqxx_no_2(v,x) - term += bielec_pxxq_no_2(x,v) * (P0tuvx_no_t(x,v,u)+P0tuvx_no_t(x,u,v)) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_inact_orb + term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & + *bielec_pxxq_no(aa,x3,v3,bb)) end do end do - term = 6.d0*term if (a.eq.b) then term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) do v=1,n_act_orb do y=1,n_act_orb do x=1,n_act_orb - term-=P0tuvx_no_t(v,x,y)*bielecCI_no(x,y,v,uu) + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt) end do end do diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f index 0f9df016..06a89318 100644 --- a/src/casscf/neworbs.irp.f +++ b/src/casscf/neworbs.irp.f @@ -24,6 +24,9 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] end do end do + do i = 1, nMonoEx + SXmatrix(i+1,i+1) += level_shift_casscf + enddo if (bavard) then do i=2,nMonoEx write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) @@ -40,73 +43,110 @@ END_PROVIDER ! Eigenvectors/eigenvalues of the single-excitation matrix END_DOC call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1) + if (bavard) then + write(6,*) ' SXdiag : lowest 5 eigenvalues ' + write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) + if(nmonoex.gt.0)then + write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2) + write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3) + write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4) + write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5) + endif + write(6,*) + write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) + endif END_PROVIDER - BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)] -&BEGIN_PROVIDER [real*8, energy_improvement] + BEGIN_PROVIDER [real*8, energy_improvement] + implicit none + if(state_following_casscf)then + energy_improvement = SXeigenval(best_vector_ovrlp_casscf) + else + energy_improvement = SXeigenval(1) + endif + END_PROVIDER + + + + BEGIN_PROVIDER [ integer, best_vector_ovrlp_casscf ] +&BEGIN_PROVIDER [ double precision, best_overlap_casscf ] + implicit none + integer :: i + double precision :: c0 + best_overlap_casscf = 0.D0 + best_vector_ovrlp_casscf = -1000 + do i=1,nMonoEx+1 + if (SXeigenval(i).lt.0.D0) then + if (abs(SXeigenvec(1,i)).gt.best_overlap_casscf) then + best_overlap_casscf=abs(SXeigenvec(1,i)) + best_vector_ovrlp_casscf = i + end if + end if + end do + if(best_vector_ovrlp_casscf.lt.0)then + best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1) + endif + c0=SXeigenvec(1,best_vector_ovrlp_casscf) + if (bavard) then + write(6,*) ' SXdiag : eigenvalue for best overlap with ' + write(6,*) ' previous orbitals = ',SXeigenval(best_vector_ovrlp_casscf) + write(6,*) ' weight of the 1st element ',c0 + endif + END_PROVIDER + + BEGIN_PROVIDER [double precision, SXvector, (nMonoEx+1)] implicit none BEGIN_DOC ! Best eigenvector of the single-excitation matrix END_DOC - integer :: ierr,matz,i - real*8 :: c0 - - if (bavard) then - write(6,*) ' SXdiag : lowest 5 eigenvalues ' - write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) - write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2) - write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3) - write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4) - write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5) - write(6,*) - write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) - endif - energy_improvement = SXeigenval(1) - - integer :: best_vector - real*8 :: best_overlap - best_overlap=0.D0 - best_vector = -1000 + integer :: i + double precision :: c0 + c0=SXeigenvec(1,best_vector_ovrlp_casscf) do i=1,nMonoEx+1 - if (SXeigenval(i).lt.0.D0) then - if (abs(SXeigenvec(1,i)).gt.best_overlap) then - best_overlap=abs(SXeigenvec(1,i)) - best_vector=i - end if - end if + SXvector(i)=SXeigenvec(i,best_vector_ovrlp_casscf)/c0 end do - if(best_vector.lt.0)then - best_vector = minloc(SXeigenval,nMonoEx+1) - endif - energy_improvement = SXeigenval(best_vector) - - c0=SXeigenvec(1,best_vector) - - if (bavard) then - write(6,*) ' SXdiag : eigenvalue for best overlap with ' - write(6,*) ' previous orbitals = ',SXeigenval(best_vector) - write(6,*) ' weight of the 1st element ',c0 - endif - - do i=1,nMonoEx+1 - SXvector(i)=SXeigenvec(i,best_vector)/c0 - end do - - -END_PROVIDER + END_PROVIDER -BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ] +BEGIN_PROVIDER [double precision, NewOrbs, (ao_num,mo_num) ] implicit none BEGIN_DOC ! Updated orbitals END_DOC integer :: i,j,ialph - call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & - NatOrbsFCI, size(NatOrbsFCI,1), & - Umat, size(Umat,1), 0.d0, & - NewOrbs, size(NewOrbs,1)) + if(state_following_casscf)then + print*,'Using the state following casscf ' + call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & + NatOrbsFCI, size(NatOrbsFCI,1), & + Umat, size(Umat,1), 0.d0, & + NewOrbs, size(NewOrbs,1)) + + level_shift_casscf *= 0.5D0 + level_shift_casscf = max(level_shift_casscf,0.002d0) + !touch level_shift_casscf + else + if(best_vector_ovrlp_casscf.ne.1.and.n_orb_swap.ne.0)then + print*,'Taking the lowest root for the CASSCF' + print*,'!!! SWAPPING MOS !!!!!!' + level_shift_casscf *= 2.D0 + level_shift_casscf = min(level_shift_casscf,0.5d0) + print*,'level_shift_casscf = ',level_shift_casscf + NewOrbs = switch_mo_coef + !mo_coef = switch_mo_coef + !soft_touch mo_coef + !call save_mos_no_occ + !stop + else + level_shift_casscf *= 0.5D0 + level_shift_casscf = max(level_shift_casscf,0.002d0) + !touch level_shift_casscf + call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & + NatOrbsFCI, size(NatOrbsFCI,1), & + Umat, size(Umat,1), 0.d0, & + NewOrbs, size(NewOrbs,1)) + endif + endif END_PROVIDER diff --git a/src/casscf/reorder_orb.irp.f b/src/casscf/reorder_orb.irp.f new file mode 100644 index 00000000..3cb90522 --- /dev/null +++ b/src/casscf/reorder_orb.irp.f @@ -0,0 +1,70 @@ +subroutine reorder_orbitals_for_casscf + implicit none + BEGIN_DOC +! routine that reorders the orbitals of the CASSCF in terms block of core, active and virtual + END_DOC + integer :: i,j,iorb + integer, allocatable :: iorder(:),array(:) + allocate(iorder(mo_num),array(mo_num)) + do i = 1, n_core_orb + iorb = list_core(i) + array(iorb) = i + enddo + + do i = 1, n_inact_orb + iorb = list_inact(i) + array(iorb) = mo_num + i + enddo + + do i = 1, n_act_orb + iorb = list_act(i) + array(iorb) = 2 * mo_num + i + enddo + + do i = 1, n_virt_orb + iorb = list_virt(i) + array(iorb) = 3 * mo_num + i + enddo + + do i = 1, mo_num + iorder(i) = i + enddo + call isort(array,iorder,mo_num) + double precision, allocatable :: mo_coef_new(:,:) + allocate(mo_coef_new(ao_num,mo_num)) + do i = 1, mo_num + mo_coef_new(:,i) = mo_coef(:,iorder(i)) + enddo + mo_coef = mo_coef_new + touch mo_coef + + list_core_reverse = 0 + do i = 1, n_core_orb + list_core(i) = i + list_core_reverse(i) = i + mo_class(i) = "Core" + enddo + + list_inact_reverse = 0 + do i = 1, n_inact_orb + list_inact(i) = i + n_core_orb + list_inact_reverse(i+n_core_orb) = i + mo_class(i+n_core_orb) = "Inactive" + enddo + + list_act_reverse = 0 + do i = 1, n_act_orb + list_act(i) = n_core_inact_orb + i + list_act_reverse(n_core_inact_orb + i) = i + mo_class(n_core_inact_orb + i) = "Active" + enddo + + list_virt_reverse = 0 + do i = 1, n_virt_orb + list_virt(i) = n_core_inact_orb + n_act_orb + i + list_virt_reverse(n_core_inact_orb + n_act_orb + i) = i + mo_class(n_core_inact_orb + n_act_orb + i) = "Virtual" + enddo + touch list_core_reverse list_core list_inact list_inact_reverse list_act list_act_reverse list_virt list_virt_reverse + +end diff --git a/src/casscf/superci_dm.irp.f b/src/casscf/superci_dm.irp.f new file mode 100644 index 00000000..0aef222b --- /dev/null +++ b/src/casscf/superci_dm.irp.f @@ -0,0 +1,207 @@ + BEGIN_PROVIDER [double precision, super_ci_dm, (mo_num,mo_num)] + implicit none + BEGIN_DOC +! density matrix of the super CI matrix, in the basis of NATURAL ORBITALS OF THE CASCI WF +! +! This is obtained from annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 +! +! WARNING ::: in the equation B3.d there is a TYPO with a forgotten MINUS SIGN (see variable mat_tmp_dm_super_ci ) + END_DOC + super_ci_dm = 0.d0 + integer :: i,j,iorb,jorb + integer :: a,aorb,b,borb + integer :: t,torb,v,vorb,u,uorb,x,xorb + double precision :: c0,ci + c0 = SXeigenvec(1,1) + ! equation B3.a of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + ! loop over the core/inact + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(iorb,iorb) = 2.d0 ! first term of B3.a + ! loop over the core/inact + do j = 1, n_core_inact_orb + jorb = list_core_inact(j) + ! loop over the virtual + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(jorb,iorb) += -2.d0 * lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,jorb) ! second term in B3.a + enddo + do t = 1, n_act_orb + torb = list_act(t) + ! thrid term of the B3.a + super_ci_dm(jorb,iorb) += - lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(jorb,torb) * (2.d0 - occ_act(t)) + enddo + enddo + enddo + + ! equation B3.b of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(iorb,torb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) + super_ci_dm(torb,iorb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(iorb,torb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + super_ci_dm(torb,iorb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + enddo + enddo + enddo + + ! equation B3.c of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(aorb,iorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) + super_ci_dm(iorb,aorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) + enddo + enddo + + ! equation B3.d of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(torb,torb) = occ_act(t) ! first term of equation B3.d + do x = 1, n_act_orb + xorb = list_act(x) + super_ci_dm(torb,torb) += - occ_act(x) * occ_act(t)* mat_tmp_dm_super_ci(x,x) ! second term involving the ONE-rdm + enddo + do u = 1, n_act_orb + uorb = list_act(u) + + ! second term of equation B3.d + do x = 1, n_act_orb + xorb = list_act(x) + do v = 1, n_act_orb + vorb = list_act(v) + super_ci_dm(torb,uorb) += 2.d0 * P0tuvx_no(v,x,t,u) * mat_tmp_dm_super_ci(v,x) ! second term involving the TWO-rdm + enddo + enddo + + ! third term of equation B3.d + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(torb,uorb) += lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(iorb,uorb) * (2.d0 - occ_act(t) - occ_act(u)) + enddo + + enddo + enddo + + ! equation B3.e of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do t = 1, n_act_orb + torb = list_act(t) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(aorb,torb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + super_ci_dm(torb,aorb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(aorb,torb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) + super_ci_dm(torb,aorb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) + enddo + enddo + enddo + + ! equation B3.f of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do a = 1, n_virt_orb + aorb = list_virt(a) + do b = 1, n_virt_orb + borb= list_virt(b) + + ! First term of equation B3.f + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(borb,aorb) += 2.d0 * lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,borb) + enddo + + ! Second term of equation B3.f + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(borb,aorb) += lowest_super_ci_coef_mo(torb,aorb) * lowest_super_ci_coef_mo(torb,borb) * occ_act(t) + enddo + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, superci_natorb, (ao_num,mo_num) +&BEGIN_PROVIDER [double precision, superci_nat_occ, (mo_num) + implicit none + call general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(super_ci_dm,mo_num,mo_num,mo_num,NatOrbsFCI,superci_nat_occ,superci_natorb) + +END_PROVIDER + + BEGIN_PROVIDER [double precision, mat_tmp_dm_super_ci, (n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC + ! computation of the term in [ ] in the equation B3.d of Roos et. al. Chemical Physics 48 (1980) 157-173 + ! + ! !!!!! WARNING !!!!!! there is a TYPO: a MINUS SIGN SHOULD APPEAR in that term + END_DOC + integer :: a,aorb,i,iorb + integer :: x,xorb,v,vorb + mat_tmp_dm_super_ci = 0.d0 + do v = 1, n_act_orb + vorb = list_act(v) + do x = 1, n_act_orb + xorb = list_act(x) + do a = 1, n_virt_orb + aorb = list_virt(a) + mat_tmp_dm_super_ci(x,v) += lowest_super_ci_coef_mo(aorb,vorb) * lowest_super_ci_coef_mo(aorb,xorb) + enddo + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + ! MARK THE MINUS SIGN HERE !!!!!!!!!!! BECAUSE OF TYPO IN THE ORIGINAL PAPER + mat_tmp_dm_super_ci(x,v) -= lowest_super_ci_coef_mo(iorb,vorb) * lowest_super_ci_coef_mo(iorb,xorb) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, lowest_super_ci_coef_mo, (mo_num,mo_num)] + implicit none + integer :: i,j,iorb,jorb + integer :: a, aorb,t, torb + double precision :: sqrt2 + + sqrt2 = 1.d0/dsqrt(2.d0) + do i = 1, nMonoEx + iorb = excit(1,i) + jorb = excit(2,i) + lowest_super_ci_coef_mo(iorb,jorb) = SXeigenvec(i+1,1) + lowest_super_ci_coef_mo(jorb,iorb) = SXeigenvec(i+1,1) + enddo + + ! a_{it} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do t = 1, n_act_orb + torb = list_act(t) + lowest_super_ci_coef_mo(torb,iorb) *= (2.d0 - occ_act(t))**(-0.5d0) + lowest_super_ci_coef_mo(iorb,torb) *= (2.d0 - occ_act(t))**(-0.5d0) + enddo + enddo + + ! a_{ia} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do a = 1, n_virt_orb + aorb = list_virt(a) + lowest_super_ci_coef_mo(aorb,iorb) *= sqrt2 + lowest_super_ci_coef_mo(iorb,aorb) *= sqrt2 + enddo + enddo + + ! a_{ta} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do a = 1, n_virt_orb + aorb = list_virt(a) + do t = 1, n_act_orb + torb = list_act(t) + lowest_super_ci_coef_mo(torb,aorb) *= occ_act(t)**(-0.5d0) + lowest_super_ci_coef_mo(aorb,torb) *= occ_act(t)**(-0.5d0) + enddo + enddo + + END_PROVIDER + diff --git a/src/casscf/swap_orb.irp.f b/src/casscf/swap_orb.irp.f new file mode 100644 index 00000000..5d442157 --- /dev/null +++ b/src/casscf/swap_orb.irp.f @@ -0,0 +1,132 @@ + BEGIN_PROVIDER [double precision, SXvector_lowest, (nMonoEx)] + implicit none + integer :: i + do i=2,nMonoEx+1 + SXvector_lowest(i-1)=SXeigenvec(i,1) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, thresh_overlap_switch] + implicit none + thresh_overlap_switch = 0.5d0 + END_PROVIDER + + BEGIN_PROVIDER [integer, max_overlap, (nMonoEx)] +&BEGIN_PROVIDER [integer, n_max_overlap] +&BEGIN_PROVIDER [integer, dim_n_max_overlap] + implicit none + double precision, allocatable :: vec_tmp(:) + integer, allocatable :: iorder(:) + allocate(vec_tmp(nMonoEx),iorder(nMonoEx)) + integer :: i + do i = 1, nMonoEx + iorder(i) = i + vec_tmp(i) = -dabs(SXvector_lowest(i)) + enddo + call dsort(vec_tmp,iorder,nMonoEx) + n_max_overlap = 0 + do i = 1, nMonoEx + if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then + n_max_overlap += 1 + max_overlap(n_max_overlap) = iorder(i) + endif + enddo + dim_n_max_overlap = max(1,n_max_overlap) + END_PROVIDER + + BEGIN_PROVIDER [integer, orb_swap, (2,dim_n_max_overlap)] +&BEGIN_PROVIDER [integer, index_orb_swap, (dim_n_max_overlap)] +&BEGIN_PROVIDER [integer, n_orb_swap ] + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer :: i,imono,iorb,jorb,j + n_orb_swap = 0 + do i = 1, n_max_overlap + imono = max_overlap(i) + iorb = excit(1,imono) + jorb = excit(2,imono) + if (excit_class(imono) == "c-a" .and.hessmat2(imono,imono).gt.0.d0)then ! core --> active rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = iorb ! core + orb_swap(2,n_orb_swap) = jorb ! active + index_orb_swap(n_orb_swap) = imono + else if (excit_class(imono) == "a-v" .and.hessmat2(imono,imono).gt.0.d0)then ! active --> virtual rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = jorb ! virtual + orb_swap(2,n_orb_swap) = iorb ! active + index_orb_swap(n_orb_swap) = imono + endif + enddo + + integer,allocatable :: orb_swap_tmp(:,:) + allocate(orb_swap_tmp(2,dim_n_max_overlap)) + do i = 1, n_orb_swap + orb_swap_tmp(1,i) = orb_swap(1,i) + orb_swap_tmp(2,i) = orb_swap(2,i) + enddo + + integer(bit_kind), allocatable :: det_i(:),det_j(:) + allocate(det_i(N_int),det_j(N_int)) + logical, allocatable :: good_orb_rot(:) + allocate(good_orb_rot(n_orb_swap)) + integer, allocatable :: index_orb_swap_tmp(:) + allocate(index_orb_swap_tmp(dim_n_max_overlap)) + index_orb_swap_tmp = index_orb_swap + good_orb_rot = .True. + integer :: icount,k + do i = 1, n_orb_swap + if(.not.good_orb_rot(i))cycle + det_i = 0_bit_kind + call set_bit_to_integer(orb_swap(1,i),det_i,N_int) + call set_bit_to_integer(orb_swap(2,i),det_i,N_int) + do j = i+1, n_orb_swap + det_j = 0_bit_kind + call set_bit_to_integer(orb_swap(1,j),det_j,N_int) + call set_bit_to_integer(orb_swap(2,j),det_j,N_int) + icount = 0 + do k = 1, N_int + icount += popcnt(ior(det_i(k),det_j(k))) + enddo + if (icount.ne.4)then + good_orb_rot(i) = .False. + good_orb_rot(j) = .False. + exit + endif + enddo + enddo + icount = n_orb_swap + n_orb_swap = 0 + do i = 1, icount + if(good_orb_rot(i))then + n_orb_swap += 1 + index_orb_swap(n_orb_swap) = index_orb_swap_tmp(i) + orb_swap(1,n_orb_swap) = orb_swap_tmp(1,i) + orb_swap(2,n_orb_swap) = orb_swap_tmp(2,i) + endif + enddo + + if(n_orb_swap.gt.0)then + print*,'n_orb_swap = ',n_orb_swap + endif + do i = 1, n_orb_swap + print*,'imono = ',index_orb_swap(i) + print*,orb_swap(1,i),'-->',orb_swap(2,i) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, switch_mo_coef, (ao_num,mo_num)] + implicit none + integer :: i,j,iorb,jorb + switch_mo_coef = NatOrbsFCI + do i = 1, n_orb_swap + iorb = orb_swap(1,i) + jorb = orb_swap(2,i) + do j = 1, ao_num + switch_mo_coef(j,jorb) = NatOrbsFCI(j,iorb) + enddo + do j = 1, ao_num + switch_mo_coef(j,iorb) = NatOrbsFCI(j,jorb) + enddo + enddo + + END_PROVIDER diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg new file mode 100644 index 00000000..5110b776 --- /dev/null +++ b/src/cipsi/EZFIO.cfg @@ -0,0 +1,5 @@ +[pert_2rdm] +type: logical +doc: If true, computes the one- and two-body rdms with perturbation theory +interface: ezfio,provider,ocaml +default: False diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 85bea747..e2917261 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -8,11 +8,6 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock] call omp_init_lock(pert_2rdm_lock) END_PROVIDER -BEGIN_PROVIDER [logical , pert_2rdm ] - implicit none - pert_2rdm = .False. -END_PROVIDER - BEGIN_PROVIDER [integer, n_orb_pert_rdm] implicit none n_orb_pert_rdm = n_act_orb diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 7825d24c..70e86fca 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -129,7 +129,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted - PROVIDE psi_det_hii N_generators_bitmask selection_weight pseudo_sym + PROVIDE psi_det_hii selection_weight pseudo_sym if (h0_type == 'SOP') then PROVIDE psi_occ_pattern_hii det_to_occ_pattern @@ -156,7 +156,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) do pt2_stoch_istate=1,N_states state_average_weight(:) = 0.d0 state_average_weight(pt2_stoch_istate) = 1.d0 - TOUCH state_average_weight pt2_stoch_istate + TOUCH state_average_weight pt2_stoch_istate selection_weight PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals pt2_w PROVIDE psi_selectors pt2_u pt2_J pt2_R @@ -523,10 +523,24 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, varianc exit else call pull_pt2_results(zmq_socket_pull, index, eI_task, vI_task, nI_task, task_id, n_tasks, b2) + if(n_tasks > pt2_n_tasks_max)then + print*,'PB !!!' + print*,'If you see this, send an email to Anthony scemama with the following content' + print*,irp_here + print*,'n_tasks,pt2_n_tasks_max = ',n_tasks,pt2_n_tasks_max + stop -1 + endif if (zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_tasks,sending) == -1) then stop 'PT2: Unable to delete tasks (send)' endif do i=1,n_tasks + if(index(i).gt.size(eI,2).or.index(i).lt.1)then + print*,'PB !!!' + print*,'If you see this, send an email to Anthony scemama with the following content' + print*,irp_here + print*,'i,index(i),size(ei,2) = ',i,index(i),size(ei,2) + stop -1 + endif eI(1:N_states, index(i)) += eI_task(1:N_states,i) vI(1:N_states, index(i)) += vI_task(1:N_states,i) nI(1:N_states, index(i)) += nI_task(1:N_states,i) @@ -706,92 +720,95 @@ END_PROVIDER - BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_W_T ] -&BEGIN_PROVIDER [ double precision, pt2_u_0 ] -&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ] - implicit none - integer :: i, t - double precision, allocatable :: tilde_w(:), tilde_cW(:) - double precision :: r, tooth_width - integer, external :: pt2_find_sample + BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ] +&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ] +&BEGIN_PROVIDER [ double precision, pt2_W_T ] +&BEGIN_PROVIDER [ double precision, pt2_u_0 ] +&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ] + implicit none + integer :: i, t + double precision, allocatable :: tilde_w(:), tilde_cW(:) + double precision :: r, tooth_width + integer, external :: pt2_find_sample + + double precision :: rss + double precision, external :: memory_of_double, memory_of_int + rss = memory_of_double(2*N_det_generators+1) + call check_mem(rss,irp_here) + + if (N_det_generators == 1) then + + pt2_w(1) = 1.d0 + pt2_cw(1) = 1.d0 + pt2_u_0 = 1.d0 + pt2_W_T = 0.d0 + pt2_n_0(1) = 0 + pt2_n_0(2) = 1 + + else + + allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) + + tilde_cW(0) = 0d0 + + do i=1,N_det_generators + tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 !+ 1.d-20 + enddo + + double precision :: norm + norm = 0.d0 + do i=N_det_generators,1,-1 + norm += tilde_w(i) + enddo + + tilde_w(:) = tilde_w(:) / norm + + tilde_cW(0) = -1.d0 + do i=1,N_det_generators + tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) + enddo + tilde_cW(:) = tilde_cW(:) + 1.d0 + + pt2_n_0(1) = 0 + do + pt2_u_0 = tilde_cW(pt2_n_0(1)) + r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth) + pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth) + if(pt2_W_T >= r - pt2_u_0) then + exit + end if + pt2_n_0(1) += 1 + if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then + print *, "teeth building failed" + stop -1 + end if + end do + + do t=2, pt2_N_teeth + r = pt2_u_0 + pt2_W_T * dble(t-1) + pt2_n_0(t) = pt2_find_sample(r, tilde_cW) + end do + pt2_n_0(pt2_N_teeth+1) = N_det_generators + + pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1)) + do t=1, pt2_N_teeth + tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) + if (tooth_width == 0.d0) then + tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + endif + ASSERT(tooth_width > 0.d0) + do i=pt2_n_0(t)+1, pt2_n_0(t+1) + pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width + end do + end do + + pt2_cW(0) = 0d0 + do i=1,N_det_generators + pt2_cW(i) = pt2_cW(i-1) + pt2_w(i) + end do + pt2_n_0(pt2_N_teeth+1) = N_det_generators - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - if (N_det_generators == 1) then - pt2_w = 1.d0 - pt2_cw = 1.d0 - pt2_W_T = 1.d0 - pt2_u_0 = 1.d0 - pt2_n_0 = 1 - return - endif - - rss = memory_of_double(2*N_det_generators+1) - call check_mem(rss,irp_here) - - allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) - - tilde_cW(0) = 0d0 - - do i=1,N_det_generators - tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 !+ 1.d-20 - enddo - - double precision :: norm - norm = 0.d0 - do i=N_det_generators,1,-1 - norm += tilde_w(i) - enddo - - tilde_w(:) = tilde_w(:) / norm - - tilde_cW(0) = -1.d0 - do i=1,N_det_generators - tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) - enddo - tilde_cW(:) = tilde_cW(:) + 1.d0 - - pt2_n_0(1) = 0 - do - pt2_u_0 = tilde_cW(pt2_n_0(1)) - r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth) - pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth) - if(pt2_W_T >= r - pt2_u_0) then - exit - end if - pt2_n_0(1) += 1 - if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then - print *, "teeth building failed" - end if - end do - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - do t=2, pt2_N_teeth - r = pt2_u_0 + pt2_W_T * dble(t-1) - pt2_n_0(t) = pt2_find_sample(r, tilde_cW) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1)) - do t=1, pt2_N_teeth - tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) - if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) - endif - ASSERT(tooth_width > 0.d0) - do i=pt2_n_0(t)+1, pt2_n_0(t+1) - pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width - end do - end do - - pt2_cW(0) = 0d0 - do i=1,N_det_generators - pt2_cW(i) = pt2_cW(i-1) + pt2_w(i) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators + endif END_PROVIDER diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 248876ef..f5217297 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -70,8 +70,6 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st) variance_match_weight(k) = product(memo_variance(k,:)) enddo - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - print *, '# var weight ', real(variance_match_weight(:),4) SOFT_TOUCH pt2_match_weight variance_match_weight end @@ -85,7 +83,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] case (0) print *, 'Using input weights in selection' - selection_weight(1:N_states) = state_average_weight(1:N_states) + selection_weight(1:N_states) = c0_weight(1:N_states) * state_average_weight(1:N_states) case (1) print *, 'Using 1/c_max^2 weight in selection' @@ -94,20 +92,30 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] case (2) print *, 'Using pt2-matching weight in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) case (3) print *, 'Using variance-matching weight in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) case (4) print *, 'Using variance- and pt2-matching weights in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states) + selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) case (5) print *, 'Using variance-matching weight in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (6) + print *, 'Using CI coefficient weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) end select + print *, '# Total weight ', real(selection_weight(:),4) END_PROVIDER @@ -165,15 +173,13 @@ subroutine select_connected(i_generator,E0,pt2,variance,norm,b,subset,csubset) call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int) - do l=1,N_generators_bitmask - do k=1,N_int - hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole,l), psi_det_generators(k,1,i_generator)) - hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator)) - particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) ) - particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) ) - enddo - call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,variance,norm,b,subset,csubset) + do k=1,N_int + hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole), psi_det_generators(k,1,i_generator)) + hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole), psi_det_generators(k,2,i_generator)) + particle_mask(k,1) = iand(generators_bitmask(k,1,s_part), not(psi_det_generators(k,1,i_generator)) ) + particle_mask(k,2) = iand(generators_bitmask(k,2,s_part), not(psi_det_generators(k,2,i_generator)) ) enddo + call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,variance,norm,b,subset,csubset) deallocate(fock_diag_tmp) end subroutine @@ -645,7 +651,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d logical :: ok integer :: s1, s2, p1, p2, ib, j, istate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef + double precision :: e_pert, delta_E, val, Hii, w, tmp, alpha_h_psi, coef double precision, external :: diag_H_mat_elem_fock double precision :: E_shift @@ -726,10 +732,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (.not.is_a_1h1p(det)) cycle endif - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - sum_e_pert = 0d0 + w = 0d0 + +! integer(bit_kind) :: occ(N_int,2), n +! call occ_pattern_of_det(det,occ,N_int) +! call occ_pattern_to_dets_size(occ,n,elec_alpha_num,N_int) do istate=1,N_states delta_E = E0(istate) - Hii + E_shift @@ -740,27 +749,43 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d tmp = -tmp endif e_pert = 0.5d0 * (tmp - delta_E) - coef = e_pert / alpha_h_psi + if (dabs(alpha_h_psi) > 1.d-4) then + coef = e_pert / alpha_h_psi + else + coef = alpha_h_psi / delta_E + endif pt2(istate) = pt2(istate) + e_pert variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi norm(istate) = norm(istate) + coef * coef - if (weight_selection /= 5) then - ! Energy selection - sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) - else - ! Variance selection - sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) - endif + + select case (weight_selection) + + case(0:4) + ! Energy selection + w = w + e_pert * selection_weight(istate) + + case(5) + ! Variance selection + w = w - alpha_h_psi * alpha_h_psi * selection_weight(istate) + + case(6) + w = w - coef * coef * selection_weight(istate) + + end select end do + + if(pseudo_sym)then - if(dabs(mat(1, p1, p2)).lt.thresh_sym)then - sum_e_pert = 10.d0 - endif + if(dabs(mat(1, p1, p2)).lt.thresh_sym)then + w = 0.d0 + endif endif - if(sum_e_pert <= buf%mini) then - call add_to_selection_buffer(buf, det, sum_e_pert) +! w = dble(n) * w + + if(w <= buf%mini) then + call add_to_selection_buffer(buf, det, w) end if end do end do diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi/selection_buffer.irp.f index 17b6e9a9..cfa3b902 100644 --- a/src/cipsi/selection_buffer.irp.f +++ b/src/cipsi/selection_buffer.irp.f @@ -198,6 +198,7 @@ subroutine make_selection_buffer_s2(b) deallocate(b%det) + print*,'n_d = ',n_d call i8sort(bit_tmp,iorder,n_d) do i=1,n_d diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 4f968ef7..b8bf6a1d 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -10,7 +10,7 @@ subroutine run_stochastic_cipsi double precision :: rss double precision, external :: memory_of_double - PROVIDE H_apply_buffer_allocated N_generators_bitmask + PROVIDE H_apply_buffer_allocated N_iter = 1 threshold_generators = 1.d0 @@ -102,7 +102,7 @@ subroutine run_stochastic_cipsi ! Add selected determinants call copy_H_apply_buffer_to_wf() - call save_wavefunction +! call save_wavefunction PROVIDE psi_coef PROVIDE psi_det diff --git a/src/cis/20.cis.bats b/src/cis/20.cis.bats index 54eefe95..bcbff701 100644 --- a/src/cis/20.cis.bats +++ b/src/cis/20.cis.bats @@ -21,6 +21,11 @@ function run() { eq $energy3 $4 $thresh } +@test "B-B" { # 2.0s + run b2_stretched.ezfio -48.995058575280950 -48.974653655601145 -48.974653655601031 + +} + @test "SiH2_3B1" { # 1.23281s 1.24958s run sih2_3b1.ezfio -289.969297318489 -289.766898643192 -289.737521023380 } diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats index 5c9ac996..d17b45a0 100644 --- a/src/cisd/30.cisd.bats +++ b/src/cisd/30.cisd.bats @@ -18,6 +18,11 @@ function run() { } +@test "B-B" { # + qp set_file b2_stretched.ezfio + run -49.120607088648597 -49.055152453388231 +} + @test "SiH2_3B1" { # 1.53842s 3.53856s qp set_file sih2_3b1.ezfio run -290.015949171697 -289.805036176618 diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f index aaa29f59..4153c9a6 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -44,8 +44,60 @@ program cisd ! * "del" orbitals which will be never occupied ! END_DOC + PROVIDE N_states read_wf = .False. SOFT_TOUCH read_wf - call only_act_bitmask - call run_cisd + call run +end + +subroutine run + implicit none + integer :: i,k + double precision :: cisdq(N_states), delta_e + double precision,external :: diag_h_mat_elem + + if(pseudo_sym)then + call H_apply_cisd_sym + else + call H_apply_cisd + endif + psi_coef = ci_eigenvectors + SOFT_TOUCH psi_coef + call save_wavefunction + call ezfio_set_cisd_energy(CI_energy) + + do i = 1,N_states + k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1) + delta_E = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int) + cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2) + enddo + print *, 'N_det = ', N_det + print*,'' + print*,'******************************' + print *, 'CISD Energies' + do i = 1,N_states + print *, i, CI_energy(i) + enddo + print*,'' + print*,'******************************' + print *, 'CISD+Q Energies' + do i = 1,N_states + print *, i, cisdq(i) + enddo + if (N_states > 1) then + print*,'' + print*,'******************************' + print*,'Excitation energies (au) (CISD+Q)' + do i = 2, N_states + print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1) + enddo + print*,'' + print*,'******************************' + print*,'Excitation energies (eV) (CISD+Q)' + do i = 2, N_states + print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, & + (cisdq(i) - cisdq(1)) / 0.0367502d0 + enddo + endif + end diff --git a/src/cisd/cisd_routine.irp.f b/src/cisd/cisd_routine.irp.f index f9e477b1..93b31e7d 100644 --- a/src/cisd/cisd_routine.irp.f +++ b/src/cisd/cisd_routine.irp.f @@ -1,17 +1,3 @@ -subroutine only_act_bitmask - implicit none - integer :: i,j,k - do k = 1, N_generators_bitmask - do j = 1, 6 - do i = 1, N_int - generators_bitmask(i,1,j,k) = act_bitmask(i,1) - generators_bitmask(i,2,j,k) = act_bitmask(i,2) - enddo - enddo - enddo - touch generators_bitmask -end - subroutine run_cisd implicit none integer :: i diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 93a91933..a8935695 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -22,6 +22,12 @@ doc: If |true|, read the wave function from the |EZFIO| file interface: ezfio,provider,ocaml default: False +[pruning] +type: float +doc: If p>0., remove p*Ndet determinants at every iteration +interface: ezfio,provider,ocaml +default: 0. + [s2_eig] type: logical doc: Force the wave function to be an eigenfunction of |S^2| @@ -32,11 +38,11 @@ default: True type: integer doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi)) interface: ezfio,provider,ocaml -default: 1 +default: 2 [weight_selection] type: integer -doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching +doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching, 6: CI coefficients interface: ezfio,provider,ocaml default: 2 diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index a930d70b..e69a1803 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -257,6 +257,18 @@ subroutine set_natural_mos double precision, allocatable :: tmp(:,:) label = "Natural" + integer :: i,j,iorb,jorb + do i = 1, n_virt_orb + iorb = list_virt(i) + do j = 1, n_core_inact_act_orb + jorb = list_core_inact_act(j) + if(one_e_dm_mo(iorb,jorb).ne. 0.d0)then + print*,'AHAHAH' + print*,iorb,jorb,one_e_dm_mo(iorb,jorb) + stop + endif + enddo + enddo call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) soft_touch mo_occ diff --git a/src/determinants/example.irp.f b/src/determinants/example.irp.f index 4d5b6b55..4f56f807 100644 --- a/src/determinants/example.irp.f +++ b/src/determinants/example.irp.f @@ -151,7 +151,7 @@ subroutine routine_example_psi_det print*,'Determinant connected' call debug_det(psi_det(1,1,idx(i)),N_int) print*,'excitation degree = ',degree_list(i) - call i_H_j(psi_det(1,1,1) , psi_det(1,1,idx(i)),hij,N_int) + call i_H_j(psi_det(1,1,1) , psi_det(1,1,idx(i)),N_int,hij) do j = 1, N_states i_H_psi(j) += hij * psi_coef(idx(i),j) enddo diff --git a/src/determinants/h_apply.irp.f b/src/determinants/h_apply.irp.f index f0d4d1c9..1c79bc75 100644 --- a/src/determinants/h_apply.irp.f +++ b/src/determinants/h_apply.irp.f @@ -124,39 +124,49 @@ subroutine copy_H_apply_buffer_to_wf PROVIDE H_apply_buffer_allocated + ASSERT (N_int > 0) ASSERT (N_det > 0) allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) ) + ! Backup determinants + j=0 do i=1,N_det - do k=1,N_int - ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) - ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num) - buffer_det(k,1,i) = psi_det(k,1,i) - buffer_det(k,2,i) = psi_det(k,2,i) - enddo + if (pruned(i)) cycle ! Pruned determinants + j+=1 + ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) + ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num) + buffer_det(:,:,j) = psi_det(:,:,i) enddo + N_det_old = j + + ! Backup coefficients do k=1,N_states + j=0 do i=1,N_det - buffer_coef(i,k) = psi_coef(i,k) + if (pruned(i)) cycle ! Pruned determinants + j += 1 + buffer_coef(j,k) = psi_coef(i,k) enddo + ASSERT ( j == N_det_old ) enddo - N_det_old = N_det + ! Update N_det + N_det = N_det_old do j=0,nproc-1 N_det = N_det + H_apply_buffer(j)%N_det enddo + ! Update array sizes if (psi_det_size < N_det) then psi_det_size = N_det TOUCH psi_det_size endif + + ! Restore backup in resized array do i=1,N_det_old - do k=1,N_int - psi_det(k,1,i) = buffer_det(k,1,i) - psi_det(k,2,i) = buffer_det(k,2,i) - enddo + psi_det(:,:,i) = buffer_det(:,:,i) ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num ) enddo @@ -165,6 +175,9 @@ subroutine copy_H_apply_buffer_to_wf psi_coef(i,k) = buffer_coef(i,k) enddo enddo + + ! Copy new buffers + !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old) & !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states,psi_det_size) diff --git a/src/determinants/h_apply_nozmq.template.f b/src/determinants/h_apply_nozmq.template.f index fac838d0..bd261bbe 100644 --- a/src/determinants/h_apply_nozmq.template.f +++ b/src/determinants/h_apply_nozmq.template.f @@ -33,22 +33,22 @@ subroutine $subroutine($params_main) do ispin=1,2 do k=1,N_int mask(k,ispin,s_hole) = & - iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,s_hole), & psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,s_part) = & - iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,s_part), & not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole1) = & - iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,d_hole1), & psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part1) = & - iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,d_part1), & not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole2) = & - iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,d_hole2), & psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part2) = & - iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + iand(generators_bitmask(k,ispin,d_part2), & not(psi_det_generators(k,ispin,i_generator)) ) enddo enddo diff --git a/src/determinants/occ_pattern.irp.f b/src/determinants/occ_pattern.irp.f index 5f37b289..6e6f9c9f 100644 --- a/src/determinants/occ_pattern.irp.f +++ b/src/determinants/occ_pattern.irp.f @@ -409,6 +409,51 @@ BEGIN_PROVIDER [ double precision, weight_occ_pattern, (N_occ_pattern,N_states) enddo END_PROVIDER +BEGIN_PROVIDER [ double precision, weight_occ_pattern_average, (N_occ_pattern) ] + implicit none + BEGIN_DOC + ! State-average weight of the occupation patterns in the wave function + END_DOC + integer :: i,j,k + weight_occ_pattern_average(:) = 0.d0 + do i=1,N_det + j = det_to_occ_pattern(i) + do k=1,N_states + weight_occ_pattern_average(j) += psi_coef(i,k) * psi_coef(i,k) * state_average_weight(k) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, psi_occ_pattern_sorted, (N_int,2,N_occ_pattern) ] +&BEGIN_PROVIDER [ double precision, weight_occ_pattern_average_sorted, (N_occ_pattern) ] +&BEGIN_PROVIDER [ integer, psi_occ_pattern_sorted_order, (N_occ_pattern) ] +&BEGIN_PROVIDER [ integer, psi_occ_pattern_sorted_order_reverse, (N_occ_pattern) ] + implicit none + BEGIN_DOC + ! Occupation patterns sorted by weight + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + allocate ( iorder(N_occ_pattern) ) + do i=1,N_occ_pattern + weight_occ_pattern_average_sorted(i) = -weight_occ_pattern_average(i) + iorder(i) = i + enddo + call dsort(weight_occ_pattern_average_sorted,iorder,N_occ_pattern) + do i=1,N_occ_pattern + do j=1,N_int + psi_occ_pattern_sorted(j,1,i) = psi_occ_pattern(j,1,iorder(i)) + psi_occ_pattern_sorted(j,2,i) = psi_occ_pattern(j,2,iorder(i)) + enddo + psi_occ_pattern_sorted_order(iorder(i)) = i + psi_occ_pattern_sorted_order_reverse(i) = iorder(i) + weight_occ_pattern_average_sorted(i) = -weight_occ_pattern_average_sorted(i) + enddo + + deallocate(iorder) + +END_PROVIDER + subroutine make_s2_eigenfunction implicit none diff --git a/src/determinants/prune_wf.irp.f b/src/determinants/prune_wf.irp.f new file mode 100644 index 00000000..c3cd8d12 --- /dev/null +++ b/src/determinants/prune_wf.irp.f @@ -0,0 +1,35 @@ +BEGIN_PROVIDER [ logical, pruned, (N_det) ] + implicit none + BEGIN_DOC + ! True if determinant is removed by pruning + END_DOC + + pruned(:) = .False. + + if (pruning == 0.d0) then + return + endif + + integer :: i,j,k,ndet_new,nsop_max + double precision :: thr + + if (s2_eig) then + + nsop_max = max(1,int ( dble(N_occ_pattern) * (1.d0 - pruning) + 0.5d0 )) + + do i=1,N_det + k = det_to_occ_pattern(i) + pruned(i) = psi_occ_pattern_sorted_order_reverse(k) > nsop_max + enddo + + else + + ndet_new = max(1,int( dble(N_det) * (1.d0 - pruning) + 0.5d0 )) + thr = psi_average_norm_contrib_sorted(ndet_new) + do i=1, N_det + pruned(i) = psi_average_norm_contrib(i) < thr + enddo + + endif + +END_PROVIDER diff --git a/src/determinants/psi_cas.irp.f b/src/determinants/psi_cas.irp.f index 8698512f..19a1c260 100644 --- a/src/determinants/psi_cas.irp.f +++ b/src/determinants/psi_cas.irp.f @@ -16,19 +16,17 @@ use bitmasks do l = 1, N_states psi_cas_coef(i,l) = 0.d0 enddo - do l=1,n_cas_bitmask - good = .True. - do k=1,N_int - good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), hf_bitmask(k,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), hf_bitmask(k,2)) ) - enddo - if (good) then - exit - endif + good = .True. + do k=1,N_int + good = good .and. ( & + iand(not(act_bitmask(k,1)), psi_det(k,1,i)) == & + iand(not(act_bitmask(k,1)), hf_bitmask(k,1)) ) .and. ( & + iand(not(act_bitmask(k,2)), psi_det(k,2,i)) == & + iand(not(act_bitmask(k,2)), hf_bitmask(k,2)) ) enddo + if (good) then + exit + endif if (good) then N_det_cas = N_det_cas+1 do k=1,N_int diff --git a/src/ezfio_files/00.create.bats b/src/ezfio_files/00.create.bats index 59bdad18..3d0eac25 100644 --- a/src/ezfio_files/00.create.bats +++ b/src/ezfio_files/00.create.bats @@ -24,6 +24,11 @@ function run { } +@test "B-B" { + qp set_file b2_stretched.ezfio + run b2_stretched.zmt 1 0 6-31g +} + @test "C2H2" { run c2h2.xyz 1 0 cc-pvdz_ecp_bfd bfd } diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index 812cd3d4..89ea3f61 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -22,7 +22,7 @@ function run_stoch() { thresh=$2 test_exe fci || skip qp set perturbation do_pt2 True - qp set determinants n_det_max 100000 + qp set determinants n_det_max $3 qp set determinants n_states 1 qp set davidson threshold_davidson 1.e-10 qp set davidson n_states_diag 1 @@ -31,137 +31,143 @@ function run_stoch() { eq $energy1 $1 $thresh } +@test "B-B" { + qp set_file b2_stretched.ezfio + qp set determinants n_det_max 10000 + qp set_frozen_core + run_stoch -49.14103054419 3.e-4 10000 +} @test "F2" { # 4.07m [[ -n $TRAVIS ]] && skip qp set_file f2.ezfio qp set_frozen_core - run_stoch -199.30486 1.e-4 + run_stoch -199.30486 1.e-4 100000 } @test "NH3" { # 10.6657s qp set_file nh3.ezfio qp set_mo_class --core="[1-4]" --act="[5-72]" - run -56.244753429144986 1.e-4 + run -56.244753429144986 1.e-4 100000 } @test "DHNO" { # 11.4721s qp set_file dhno.ezfio qp set_mo_class --core="[1-7]" --act="[8-64]" - run -130.459020029816 1.e-4 + run -130.459020029816 1.e-4 100000 } @test "HCO" { # 12.2868s qp set_file hco.ezfio - run -113.297494345682 1.e-4 + run -113.297494345682 1.e-4 100000 } @test "H2O2" { # 12.9214s qp set_file h2o2.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]" - run -151.00477 1.e-4 + run -151.00477 1.e-4 100000 } @test "HBO" { # 13.3144s [[ -n $TRAVIS ]] && skip qp set_file hbo.ezfio - run -100.212829869715 1.e-4 + run -100.212829869715 1.e-4 100000 } @test "H2O" { # 11.3727s [[ -n $TRAVIS ]] && skip qp set_file h2o.ezfio - run -76.2359268957699 1.e-4 + run -76.2359268957699 1.e-4 100000 } @test "ClO" { # 13.3755s [[ -n $TRAVIS ]] && skip qp set_file clo.ezfio - run -534.545881614967 1.e-4 + run -534.545881614967 1.e-4 100000 } @test "SO" { # 13.4952s [[ -n $TRAVIS ]] && skip qp set_file so.ezfio - run -26.0158153138924 1.e-4 + run -26.0126927641744 1.e-4 100000 } @test "H2S" { # 13.6745s [[ -n $TRAVIS ]] && skip qp set_file h2s.ezfio - run -398.859168655255 1.e-4 + run -398.859168655255 1.e-4 100000 } @test "OH" { # 13.865s [[ -n $TRAVIS ]] && skip qp set_file oh.ezfio - run -75.6120779012574 1.e-4 + run -75.6120779012574 1.e-4 100000 } @test "SiH2_3B1" { # 13.938ss [[ -n $TRAVIS ]] && skip qp set_file sih2_3b1.ezfio - run -290.017539006762 1.e-4 + run -290.017539006762 1.e-4 100000 } @test "H3COH" { # 14.7299s [[ -n $TRAVIS ]] && skip qp set_file h3coh.ezfio - run -115.205941463667 1.e-4 + run -115.205941463667 1.e-4 100000 } @test "SiH3" { # 15.99s [[ -n $TRAVIS ]] && skip qp set_file sih3.ezfio - run -5.57241217753818 1.e-4 + run -5.57241217753818 1.e-4 100000 } @test "CH4" { # 16.1612s [[ -n $TRAVIS ]] && skip qp set_file ch4.ezfio qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" - run -40.2409678239136 1.e-4 + run -40.2409678239136 1.e-4 100000 } @test "ClF" { # 16.8864s [[ -n $TRAVIS ]] && skip qp set_file clf.ezfio - run -559.170272077166 1.e-4 + run -559.170128224959 1.e-4 100000 } @test "SO2" { # 17.5645s [[ -n $TRAVIS ]] && skip qp set_file so2.ezfio qp set_mo_class --core="[1-8]" --act="[9-87]" - run -41.5746738713298 1.e-4 + run -41.5746738713298 1.e-4 100000 } @test "C2H2" { # 17.6827s [[ -n $TRAVIS ]] && skip qp set_file c2h2.ezfio qp set_mo_class --act="[1-30]" --del="[31-36]" - run -12.3656179738175 1.e-4 + run -12.3658547549095 1.e-4 100000 } @test "N2" { # 18.0198s [[ -n $TRAVIS ]] && skip qp set_file n2.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]" - run -109.291600196629 1.e-4 + run -109.291711886659 1.e-4 100000 } @test "N2H4" { # 18.5006s [[ -n $TRAVIS ]] && skip qp set_file n2h4.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]" - run -111.367332681559 1.e-4 + run -111.367332681559 1.e-4 100000 } @test "CO2" { # 21.1748s [[ -n $TRAVIS ]] && skip qp set_file co2.ezfio qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]" - run -187.968599504402 1.e-4 + run -187.968599504402 1.e-4 100000 } @@ -169,13 +175,13 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file cu_nh3_4_2plus.ezfio qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]" - run -1862.98614665139 1.e-04 + run -1862.98614665139 1.e-04 100000 } @test "HCN" { # 20.3273s [[ -n $TRAVIS ]] && skip qp set_file hcn.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]" - run -93.0728641601823 1.e-4 + run -93.0791660745576 1.e-4 100000 } diff --git a/src/generators_fluid/NEED b/src/generators_fluid/NEED deleted file mode 100644 index d3d4d2c7..00000000 --- a/src/generators_fluid/NEED +++ /dev/null @@ -1 +0,0 @@ -determinants diff --git a/src/generators_fluid/README.rst b/src/generators_fluid/README.rst deleted file mode 100644 index e69de29b..00000000 diff --git a/src/generators_fluid/extract_cas.irp.f b/src/generators_fluid/extract_cas.irp.f deleted file mode 100644 index 9cdaf27f..00000000 --- a/src/generators_fluid/extract_cas.irp.f +++ /dev/null @@ -1,23 +0,0 @@ -subroutine extract_cas - implicit none - BEGIN_DOC - ! Replaces the total wave function by the normalized projection on the CAS. - END_DOC - - integer :: i,j,k - do k=1,N_states - do j=1,N_det_generators - psi_coef(j,k) = psi_coef_generators(j,k) - enddo - enddo - - do j=1,N_det_generators - do k=1,N_int - psi_det(k,1,j) = psi_det_generators(k,1,j) - psi_det(k,2,j) = psi_det_generators(k,2,j) - enddo - enddo - N_det = N_det_generators - - SOFT_TOUCH N_det psi_det psi_coef -end diff --git a/src/generators_fluid/generators.irp.f b/src/generators_fluid/generators.irp.f deleted file mode 100644 index 153ab605..00000000 --- a/src/generators_fluid/generators.irp.f +++ /dev/null @@ -1,101 +0,0 @@ -use bitmasks - -BEGIN_PROVIDER [ character*(32), generators_type] - implicit none - generators_type = trim("CAS") - -END_PROVIDER - -BEGIN_PROVIDER [ integer, N_det_generators ] - implicit none - BEGIN_DOC - ! Number of generator detetrminants - END_DOC - if(generators_type == "CAS")then - N_det_generators = N_det_generators_CAS - else if (generators_type == "HF")then - N_det_generators = N_det_generators_HF - else if (generators_type == "HF_SD")then - N_det_generators = N_det_generators_HF_SD - endif - N_det_generators = max(N_det_generators,1) - call write_int(6,N_det_generators,'Number of generators') -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - - if(generators_type == "CAS")then - psi_det_generators(1:N_int,1:2,1:N_det_generators_CAS) = psi_det_generators_CAS(1:N_int,1:2,1:N_det_generators_CAS) - psi_coef_generators(1:N_det_generators_CAS,1:N_states) = psi_coef_generators_CAS(1:N_det_generators_CAS,1:N_states) - else if (generators_type == "HF")then - psi_det_generators(1:N_int,1:2,1:N_det_generators_HF) = psi_det_generators_HF(1:N_int,1:2,1:N_det_generators_HF) - psi_coef_generators(1:N_det_generators_HF,1:N_states) = psi_coef_generators_HF(1:N_det_generators_HF,1:N_states) - else if (generators_type == "HF_SD")then - psi_det_generators(1:N_int,1:2,1:N_det_generators_HF_SD) = psi_det_generators_HF_SD(1:N_int,1:2,1:N_det_generators_HF_SD) - psi_coef_generators(1:N_det_generators_HF_SD,1:N_states) = psi_coef_generators_HF_SD(1:N_det_generators_HF_SD,1:N_states) - endif - -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order, (psi_det_size) ] - - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - if(generators_type == "CAS")then - psi_det_sorted_gen = psi_det_sorted_gen_CAS - psi_coef_sorted_gen = psi_coef_sorted_gen_CAS - psi_det_sorted_gen_order = psi_det_sorted_gen_CAS_order - else if(generators_type == "HF")then - psi_det_sorted_gen = 0_bit_kind - psi_coef_sorted_gen = 0.d0 - psi_det_sorted_gen_order = 0 - else if(generators_type == "HF_SD")then - psi_det_sorted_gen = psi_det_sorted_gen_HF_SD - psi_coef_sorted_gen = psi_coef_sorted_gen_HF_SD - psi_det_sorted_gen_order = psi_det_sorted_gen_HF_SD_order - endif -END_PROVIDER - - -BEGIN_PROVIDER [integer, degree_max_generators] - implicit none - BEGIN_DOC -! Max degree of excitation (respect to HF) of the generators - END_DOC - integer :: i,degree - degree_max_generators = 0 - do i = 1, N_det_generators - call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int) - if(degree .gt. degree_max_generators)then - degree_max_generators = degree - endif - enddo -END_PROVIDER - -BEGIN_PROVIDER [ integer, size_select_max] - implicit none - BEGIN_DOC - ! Size of the select_max array - END_DOC - size_select_max = 10000 -END_PROVIDER - -BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] - implicit none - BEGIN_DOC - ! Memo to skip useless selectors - END_DOC - select_max = huge(1.d0) -END_PROVIDER - diff --git a/src/generators_fluid/generators_cas.irp.f b/src/generators_fluid/generators_cas.irp.f deleted file mode 100644 index b6d83e0a..00000000 --- a/src/generators_fluid/generators_cas.irp.f +++ /dev/null @@ -1,69 +0,0 @@ -use bitmasks - -BEGIN_PROVIDER [ integer, N_det_generators_CAS ] - implicit none - BEGIN_DOC - ! Number of generator detetrminants - END_DOC - integer :: i,k,l - logical :: good - integer, external :: number_of_holes,number_of_particles - call write_time(6) - N_det_generators_CAS = 0 - do i=1,N_det - good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 ) - if (good) then - N_det_generators_CAS += 1 - endif - enddo - N_det_generators_CAS = max(N_det_generators_CAS,1) - call write_int(6,N_det_generators_CAS,'Number of generators_CAS') -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_CAS, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators_CAS, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen_CAS, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen_CAS, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_CAS_order, (psi_det_size) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the gen_CASerator is the - ! Hartree-Fock determinant - END_DOC - integer :: i, k, l, m - logical :: good - integer, external :: number_of_holes,number_of_particles - integer, allocatable :: nongen_CAS(:) - integer :: inongen_CAS - - allocate(nongen_CAS(N_det)) - - inongen_CAS = 0 - m=0 - do i=1,N_det - good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 ) - if (good) then - m = m+1 - psi_det_sorted_gen_CAS_order(i) = m - do k=1,N_int - psi_det_generators_CAS(k,1,m) = psi_det_sorted(k,1,i) - psi_det_generators_CAS(k,2,m) = psi_det_sorted(k,2,i) - enddo - psi_coef_generators_CAS(m,:) = psi_coef_sorted(i,:) - else - inongen_CAS += 1 - nongen_CAS(inongen_CAS) = i - endif - enddo - ASSERT (m == N_det_generators_CAS) - - psi_det_sorted_gen_CAS(:,:,:N_det_generators_CAS) = psi_det_generators_CAS(:,:,:N_det_generators_CAS) - psi_coef_sorted_gen_CAS(:N_det_generators_CAS, :) = psi_coef_generators_CAS(:N_det_generators_CAS, :) - do i=1,inongen_CAS - psi_det_sorted_gen_CAS_order(nongen_CAS(i)) = N_det_generators_CAS+i - psi_det_sorted_gen_CAS(:,:,N_det_generators_CAS+i) = psi_det_sorted(:,:,nongen_CAS(i)) - psi_coef_sorted_gen_CAS(N_det_generators_CAS+i, :) = psi_coef_sorted(nongen_CAS(i),:) - end do - -END_PROVIDER - diff --git a/src/generators_fluid/generators_hf.irp.f b/src/generators_fluid/generators_hf.irp.f deleted file mode 100644 index d4d2e728..00000000 --- a/src/generators_fluid/generators_hf.irp.f +++ /dev/null @@ -1,51 +0,0 @@ - -use bitmasks - -BEGIN_PROVIDER [ integer, N_det_generators_HF ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the number of generators is 1 : the - ! Hartree-Fock determinant - END_DOC - N_det_generators_HF = 1 -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_HF, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators_HF, (psi_det_size,N_states) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - psi_det_generators_HF = 0_bit_kind - integer :: i,j - integer :: degree - - do i=1,N_int - psi_det_generators_HF(i,1,1) = HF_bitmask(i,1) - psi_det_generators_HF(i,2,1) = HF_bitmask(i,2) - enddo - - do j=1,N_det - call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int) - if (degree == 0) then - exit - endif - end do - - psi_det_generators_HF(:,:,1) = psi_det(:,:,j) - psi_coef_generators_HF(1,:) = 1.d0 - -END_PROVIDER - - BEGIN_PROVIDER [ integer , HF_index ] - implicit none - integer :: j,degree - do j=1,N_det - call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,j),degree,N_int) - if (degree == 0) then - HF_index = j - exit - endif - end do -END_PROVIDER diff --git a/src/generators_fluid/generators_hf_sd.irp.f b/src/generators_fluid/generators_hf_sd.irp.f deleted file mode 100644 index 9c13a5a0..00000000 --- a/src/generators_fluid/generators_hf_sd.irp.f +++ /dev/null @@ -1,80 +0,0 @@ - -use bitmasks - -BEGIN_PROVIDER [ integer, N_det_generators_HF_SD ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the number of generators is 1 : the - ! Hartree-Fock determinant - END_DOC - N_det_generators_HF_SD = 0 - integer :: i,degree - double precision :: thr - double precision :: accu - accu = 0.d0 - thr = threshold_generators - do i = 1, N_det - call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,i),degree,N_int) - if(degree.le.2.and. accu .le. thr )then - accu += psi_coef_sorted(i,1)**2 - N_det_generators_HF_SD += 1 - endif - enddo -!print*,'' -!print*,'N_det_generators_HF_SD = ',N_det_generators_HF_SD -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_HF_SD, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators_HF_SD, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen_HF_SD, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen_HF_SD, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_HF_SD_order, (psi_det_size) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - psi_det_generators_HF_SD = 0_bit_kind - integer :: i,j,k - integer :: degree - double precision :: thr - double precision :: accu - integer, allocatable :: nongen(:) - integer :: inongen - - allocate(nongen(N_det)) - - thr = threshold_generators - - accu = 0.d0 - k = 0 - inongen = 0 - do j=1,N_det - call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,j),degree,N_int) - if(degree.le.2.and. accu.le.thr )then - accu += psi_coef_sorted(j,1)**2 - k += 1 - psi_det_sorted_gen_HF_SD_order(j) = k - do i = 1, N_int - psi_det_generators_HF_SD(i,1,k) = psi_det_sorted(i,1,j) - psi_det_generators_HF_SD(i,2,k) = psi_det_sorted(i,2,j) - enddo - do i = 1, N_states - psi_coef_generators_HF_SD(k,i) = psi_coef_sorted(j,i) - enddo - else - inongen += 1 - nongen(inongen) = j - endif - end do - - psi_det_sorted_gen_HF_SD(:,:,:N_det_generators_HF_SD) = psi_det_generators_HF_SD(:,:,:N_det_generators_HF_SD) - psi_coef_sorted_gen_HF_SD(:N_det_generators_HF_SD, :) = psi_coef_generators_HF_SD(:N_det_generators_HF_SD, :) - do i=1,inongen - psi_det_sorted_gen_HF_SD_order(nongen(i)) = N_det_generators_HF_SD+i - psi_det_sorted_gen_HF_SD(:,:,N_det_generators_HF_SD+i) = psi_det_sorted(:,:,nongen(i)) - psi_coef_sorted_gen_HF_SD(N_det_generators_HF_SD+i, :) = psi_coef_sorted(nongen(i),:) - end do - -END_PROVIDER - diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index ae78309a..8a9dde37 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -17,6 +17,10 @@ function run() { } +@test "B-B" { # 3s + run b2_stretched.ezfio -48.9950585752809 +} + @test "SiH2_3B1" { # 0.539000 1.51094s run sih2_3b1.ezfio -289.9654718650881 } diff --git a/src/kohn_sham_rs/61.rsks.bats b/src/kohn_sham_rs/61.rsks.bats index 558c5027..c5e67350 100644 --- a/src/kohn_sham_rs/61.rsks.bats +++ b/src/kohn_sham_rs/61.rsks.bats @@ -21,7 +21,6 @@ function run() { eq $energy $3 $thresh } - @test "H3COH" { run h3coh.ezfio sr_pbe -115.50238225208 } diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index 126705bf..a055aad3 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -23,7 +23,7 @@ size: (mo_basis.mo_num) [mo_class] type: MO_class doc: [ Core | Inactive | Active | Virtual | Deleted ], as defined by :ref:`qp_set_mo_class` -interface: ezfio, provider +interface: ezfio size: (mo_basis.mo_num) [ao_md5] diff --git a/src/mo_basis/mo_class.irp.f b/src/mo_basis/mo_class.irp.f new file mode 100644 index 00000000..95fbb443 --- /dev/null +++ b/src/mo_basis/mo_class.irp.f @@ -0,0 +1,40 @@ +! DO NOT MODIFY BY HAND +! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py +! from file /home/eginer/programs/qp2/src/mo_basis/EZFIO.cfg + + +BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ] + implicit none + BEGIN_DOC +! [ Core | Inactive | Active | Virtual | Deleted ], as defined by :ref:`qp_set_mo_class` + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + if (size(mo_class) == 0) return + + call ezfio_has_mo_basis_mo_class(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: mo_class ] <<<<< ..' + call ezfio_get_mo_basis_mo_class(mo_class) + else + mo_class(:) = 'Active' + endif + 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( mo_class, (mo_num)*32, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_class with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 610e9a8c..aa04fb01 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -91,7 +91,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ] enddo enddo endif - END_PROVIDER BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ] diff --git a/src/mo_basis/utils.irp.f b/src/mo_basis/utils.irp.f index e141867a..12c6c79d 100644 --- a/src/mo_basis/utils.irp.f +++ b/src/mo_basis/utils.irp.f @@ -4,7 +4,6 @@ subroutine save_mos integer :: i,j call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) - call ezfio_set_mo_basis_mo_num(mo_num) call ezfio_set_mo_basis_mo_label(mo_label) call ezfio_set_mo_basis_ao_md5(ao_md5) @@ -17,6 +16,29 @@ subroutine save_mos enddo call ezfio_set_mo_basis_mo_coef(buffer) call ezfio_set_mo_basis_mo_occ(mo_occ) + call ezfio_set_mo_basis_mo_class(mo_class) + deallocate (buffer) + +end + + +subroutine save_mos_no_occ + implicit none + double precision, allocatable :: buffer(:,:) + integer :: i,j + + call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) + !call ezfio_set_mo_basis_mo_num(mo_num) + !call ezfio_set_mo_basis_mo_label(mo_label) + !call ezfio_set_mo_basis_ao_md5(ao_md5) + allocate ( buffer(ao_num,mo_num) ) + buffer = 0.d0 + do j = 1, mo_num + do i = 1, ao_num + buffer(i,j) = mo_coef(i,j) + enddo + enddo + call ezfio_set_mo_basis_mo_coef(buffer) deallocate (buffer) end @@ -40,6 +62,7 @@ subroutine save_mos_truncated(n) enddo call ezfio_set_mo_basis_mo_coef(buffer) call ezfio_set_mo_basis_mo_occ(mo_occ) + call ezfio_set_mo_basis_mo_class(mo_class) deallocate (buffer) end @@ -217,3 +240,64 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label) end +subroutine mo_coef_new_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,mo_coef_before,eig,mo_coef_new) + implicit none + BEGIN_DOC +! You enter with matrix in the MO basis defined with the mo_coef_before. +! +! You SVD the matrix and set the eigenvectors as mo_coef_new ordered by increasing singular values + END_DOC + integer,intent(in) :: lda,m,n + double precision, intent(in) :: matrix(lda,n),mo_coef_before(ao_num,m) + double precision, intent(out) :: eig(m),mo_coef_new(ao_num,m) + + integer :: i,j + double precision :: accu + double precision, allocatable :: mo_coef_tmp(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, Vt, A + + call write_time(6) + if (m /= mo_num) then + print *, irp_here, ': Error : m/= mo_num' + stop 1 + endif + + allocate(A(lda,n),U(lda,n),D(m),Vt(lda,n),mo_coef_tmp(ao_num,mo_num)) + + do j=1,n + do i=1,m + A(i,j) = matrix(i,j) + enddo + enddo + mo_coef_tmp = mo_coef_before + + call svd(A,lda,U,lda,D,Vt,lda,m,n) + + write (6,'(A)') '' + write (6,'(A)') 'Eigenvalues' + write (6,'(A)') '-----------' + write (6,'(A)') '' + write (6,'(A)') '======== ================ ================' + write (6,'(A)') ' MO Eigenvalue Cumulative ' + write (6,'(A)') '======== ================ ================' + + accu = 0.d0 + do i=1,m + accu = accu + D(i) + write (6,'(I8,1X,F16.10,1X,F16.10)') i,D(i), accu + enddo + write (6,'(A)') '======== ================ ================' + write (6,'(A)') '' + + call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_tmp,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef_new,size(mo_coef_new,1)) + + do i=1,m + eig(i) = D(i) + enddo + + deallocate(A,U,Vt,D,mo_coef_tmp) + call write_time(6) + +end + + diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index 6bbbfa39..f70ed0de 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -6,6 +6,7 @@ program molden character*(128) :: output integer :: i_unit_output,getUnitAndOpen integer :: i,j,k,l + double precision, parameter :: a0 = 0.529177249d0 PROVIDE ezfio_filename @@ -22,7 +23,7 @@ program molden trim(element_name(int(nucl_charge(i)))), & i, & int(nucl_charge(i)), & - nucl_coord(i,1), nucl_coord(i,2), nucl_coord(i,3) + nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0 enddo write(i_unit_output,'(A)') '[GTO]' diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f index 3323b46e..a92d1a51 100644 --- a/src/tools/print_wf.irp.f +++ b/src/tools/print_wf.irp.f @@ -14,7 +14,7 @@ program print_wf ! this has to be done in order to be sure that N_det, psi_det and - ! psi_coef are the wave function stored in the |EZFIO| directory. + ! psi_coef_sorted are the wave function stored in the |EZFIO| directory. read_wf = .True. touch read_wf call routine @@ -45,15 +45,15 @@ subroutine routine do i = 1, min(N_det_print_wf,N_det) print*,'' print*,'i = ',i - call debug_det(psi_det(1,1,i),N_int) - call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,1),degree,N_int) + call debug_det(psi_det_sorted(1,1,i),N_int) + call get_excitation_degree(psi_det_sorted(1,1,i),psi_det_sorted(1,1,1),degree,N_int) print*,'degree = ',degree if(degree == 0)then print*,'Reference determinant ' - call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,h00) + call i_H_j(psi_det_sorted(1,1,i),psi_det_sorted(1,1,i),N_int,h00) else if(degree .le. 2)then - call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii) - call i_H_j(psi_det(1,1,1),psi_det(1,1,i),N_int,hij) + call i_H_j(psi_det_sorted(1,1,i),psi_det_sorted(1,1,i),N_int,hii) + call i_H_j(psi_det_sorted(1,1,1),psi_det_sorted(1,1,i),N_int,hij) delta_e = hii - h00 coef_1 = hij/(h00-hii) if(hij.ne.0.d0)then @@ -65,25 +65,25 @@ subroutine routine else coef_2_2 = 0.d0 endif - call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int) + call get_excitation(psi_det_sorted(1,1,1),psi_det_sorted(1,1,i),exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) print*,'phase = ',phase if(degree == 1)then print*,'s1',s1 print*,'h1,p1 = ',h1,p1 if(s1 == 1)then - norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1)) - norm_mono_a_2 += dabs(psi_coef(i,1)/psi_coef(1,1))**2 + norm_mono_a += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1)) + norm_mono_a_2 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))**2 norm_mono_a_pert += dabs(coef_1) norm_mono_a_pert_2 += dabs(coef_1)**2 else - norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1)) - norm_mono_b_2 += dabs(psi_coef(i,1)/psi_coef(1,1))**2 + norm_mono_b += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1)) + norm_mono_b_2 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))**2 norm_mono_b_pert += dabs(coef_1) norm_mono_b_pert_2 += dabs(coef_1)**2 endif double precision :: hmono,hdouble - call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble,phase) + call i_H_j_verbose(psi_det_sorted(1,1,1),psi_det_sorted(1,1,i),N_int,hij,hmono,hdouble,phase) print*,'hmono = ',hmono print*,'hdouble = ',hdouble print*,'hmono+hdouble = ',hmono+hdouble @@ -99,9 +99,9 @@ subroutine routine print*,'Delta E = ',h00-hii print*,'coef pert (1) = ',coef_1 print*,'coef 2x2 = ',coef_2_2 - print*,'Delta E_corr = ',psi_coef(i,1)/psi_coef(1,1) * hij + print*,'Delta E_corr = ',psi_coef_sorted(i,1)/psi_coef_sorted(1,1) * hij endif - print*,'amplitude = ',psi_coef(i,1)/psi_coef(1,1) + print*,'amplitude = ',psi_coef_sorted(i,1)/psi_coef_sorted(1,1) enddo diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index d441e1df..2bcd04dc 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -9,7 +9,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 1 @@ -26,7 +26,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 2 @@ -43,7 +43,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin print*,'' @@ -70,7 +70,7 @@ END_DOC double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 4 @@ -79,7 +79,9 @@ double precision :: wall_0,wall_1 call wall_time(wall_0) print*,'providing the state average TWO-RDM ...' - call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + print*,'psi_det_size = ',psi_det_size + print*,'N_det = ',N_det + call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,N_states,size(psi_coef,1)) call wall_time(wall_1) print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 diff --git a/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f b/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f index 386e2a54..baa26ced 100644 --- a/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f @@ -7,7 +7,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 1 @@ -24,7 +24,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 2 @@ -41,7 +41,7 @@ ! = END_DOC allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin print*,'' @@ -68,7 +68,7 @@ END_DOC double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) + state_weights = state_average_weight integer :: ispin ! condition for alpha/beta spin ispin = 4 diff --git a/src/two_body_rdm/orb_range_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f index a8684185..058ed1c5 100644 --- a/src/two_body_rdm/orb_range_routines.irp.f +++ b/src/two_body_rdm/orb_range_routines.irp.f @@ -147,6 +147,7 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l do i=1,maxab idx0(i) = i enddo + ! Prepare the array of all alpha single excitations ! ------------------------------------------------- diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index ed9932c9..2a655eed 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -57,6 +57,8 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j + + if(isize<2)return c = x( shiftr(first+last,1) ) i = first diff --git a/tests/input/b2_stretched.zmt b/tests/input/b2_stretched.zmt new file mode 100644 index 00000000..04950a9b --- /dev/null +++ b/tests/input/b2_stretched.zmt @@ -0,0 +1,3 @@ +b +b 1 3.0 + diff --git a/tests/input/n2_stretched.xyz b/tests/input/n2_stretched.xyz new file mode 100644 index 00000000..28e26041 --- /dev/null +++ b/tests/input/n2_stretched.xyz @@ -0,0 +1,4 @@ +2 +N2 stretched +N 0. 0. 0. +N 0. 0. 2.1167090