diff --git a/.drone.yml b/.drone.yml new file mode 100644 index 00000000..d3d3ef92 --- /dev/null +++ b/.drone.yml @@ -0,0 +1,52 @@ +--- +kind: pipeline +type: docker +name: gfortran-debug + +clone: + depth: 10 + +steps: + - name: configure debug + image: scemama666/qp2_env + commands: + - ./configure -i all -c ./config/gfortran_debug.cfg + - bash -c "source quantum_package.rc ; exec qp_plugins download https://gitlab.com/scemama/qp_plugins_scemama" + - bash -c "source quantum_package.rc ; exec qp_plugins install champ" + + - name: compile debug + image: scemama666/qp2_env + commands: + - bash -c "source quantum_package.rc ; exec ninja" + + - name: testing debug + image: scemama666/qp2_env + commands: + - bash -c "source quantum_package.rc ; TRAVIS=1 exec qp_test -a" + + - name: configure fast + image: scemama666/qp2_env + commands: + - ./configure -c ./config/gfortran_avx.cfg + + - name: compile fast + image: scemama666/qp2_env + commands: + - bash -c "source quantum_package.rc ; exec ninja" + + - name: testing fast + image: scemama666/qp2_env + commands: + - bash -c "source quantum_package.rc ; exec qp_test -a" + + - name: notify + image: drillster/drone-email + settings: + host: + from_secret: hostname # irsamc.ups-tlse.fr + from: + from_secret: from # drone@irssv7.ups-tlse.fr + recipients: + from_secret: recipients # scemama@irsamc.ups-tlse.fr + when: + status: [changed, failure] diff --git a/README.md b/README.md index b8141b88..24d6277e 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # Quantum Package 2.2 - + + [![DOI](https://zenodo.org/badge/167513335.svg)](https://zenodo.org/badge/latestdoi/167513335) @@ -44,9 +45,19 @@ https://arxiv.org/abs/1902.08154 # Credits +* [TREX Center of Excellence](https://trex-coe.eu) +* [ERC PTEROSOR](https://lcpq.github.io/PTEROSOR) * [CNRS](http://www.cnrs.fr) * [Laboratoire de Chimie et Physique Quantiques](http://lcpq.ups-tlse.fr) * [Laboratoire de Chimie Théorique](http://www.lct.jussieu.fr) * [Argonne Leadership Computing Facility](http://alcf.anl.gov) * [CALMIP](https://www.calmip.univ-toulouse.fr) + +------------------------------ + + + +[TREX: Targeting Real Chemical Accuracy at the Exascale](https://trex-coe.eu) project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content. + +[PTEROSOR](https://lcpq.github.io/PTEROSOR) project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863481). diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 091423e4..e7c44b37 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -146,6 +146,17 @@ def write_ezfio(res, filename): ezfio.set_ao_basis_ao_nucl(at) ezfio.set_ao_basis_ao_prim_num(num_prim) ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) + try: + normf = res.normf + if normf == 0: + ezfio.set_ao_basis_ao_normalized(True) + elif normf == 1: + ezfio.set_ao_basis_ao_normalized(False) + else: + print("BUG in NORMF") + sys.exit(0) + except AttributeError: + ezfio.set_ao_basis_ao_normalized(True) # ~#~#~#~#~#~#~ # # P a r s i n g # @@ -224,7 +235,7 @@ def write_ezfio(res, filename): exponent += [p.expo for p in b.prim] ang_mom.append(str.count(s, "z")) shell_prim_num.append(len(b.prim)) - shell_index += [nshell_tot+1] * len(b.prim) + shell_index += [nshell_tot] * len(b.prim) # ~#~#~#~#~ # # W r i t e # diff --git a/bin/qp_test b/bin/qp_test index 67b3ea02..288b7291 100755 --- a/bin/qp_test +++ b/bin/qp_test @@ -60,19 +60,14 @@ def main(arguments): print("Running tests for %s"%(bats_file)) print("") if arguments["-v"]: - p = None if arguments["TEST"]: test = "export TEST=%s ; "%arguments["TEST"] else: test = "" - try: - os.system(test+" python3 bats_to_sh.py "+bats_file+ + os.system(test+" python3 bats_to_sh.py "+bats_file+ "| bash") - except: - if p: - p.terminate() else: - subprocess.check_call(["bats", bats_file], env=os.environ) + subprocess.check_call(["bats", "--verbose-run", "--trace", bats_file], env=os.environ) diff --git a/config/ifort_2019_avx.cfg b/config/ifort_2019_avx.cfg index 661a0e8f..c5bed0d8 100644 --- a/config/ifort_2019_avx.cfg +++ b/config/ifort_2019_avx.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_avx_mpi.cfg b/config/ifort_2019_avx_mpi.cfg index 2d212db5..5b4d2922 100644 --- a/config/ifort_2019_avx_mpi.cfg +++ b/config/ifort_2019_avx_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_mpi_rome.cfg b/config/ifort_2019_mpi_rome.cfg index 171219e6..054d4d7d 100644 --- a/config/ifort_2019_mpi_rome.cfg +++ b/config/ifort_2019_mpi_rome.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_rome.cfg b/config/ifort_2019_rome.cfg index e923a1dd..a18a0acb 100644 --- a/config/ifort_2019_rome.cfg +++ b/config/ifort_2019_rome.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_sse4.cfg b/config/ifort_2019_sse4.cfg index a3aa7cbd..2cdbc2c5 100644 --- a/config/ifort_2019_sse4.cfg +++ b/config/ifort_2019_sse4.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_sse4_mpi.cfg b/config/ifort_2019_sse4_mpi.cfg index 6959d176..d20cd2a2 100644 --- a/config/ifort_2019_sse4_mpi.cfg +++ b/config/ifort_2019_sse4_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED diff --git a/config/ifort_2019_xHost.cfg b/config/ifort_2019_xHost.cfg index 22d28803..59c6146b 100644 --- a/config/ifort_2019_xHost.cfg +++ b/config/ifort_2019_xHost.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=64 -DINTEL -DSET_NESTED diff --git a/config/ifort_2021_avx.cfg b/config/ifort_2021_avx.cfg index 6f657052..6c34cf47 100644 --- a/config/ifort_2021_avx.cfg +++ b/config/ifort_2021_avx.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_avx_mpi.cfg b/config/ifort_2021_avx_mpi.cfg index c991a4a9..4c893c73 100644 --- a/config/ifort_2021_avx_mpi.cfg +++ b/config/ifort_2021_avx_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_mpi_rome.cfg b/config/ifort_2021_mpi_rome.cfg index 8413d23d..e47a466e 100644 --- a/config/ifort_2021_mpi_rome.cfg +++ b/config/ifort_2021_mpi_rome.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_rome.cfg b/config/ifort_2021_rome.cfg index b3023186..504438c9 100644 --- a/config/ifort_2021_rome.cfg +++ b/config/ifort_2021_rome.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_sse4.cfg b/config/ifort_2021_sse4.cfg index a6299665..07c3ebb8 100644 --- a/config/ifort_2021_sse4.cfg +++ b/config/ifort_2021_sse4.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_sse4_mpi.cfg b/config/ifort_2021_sse4_mpi.cfg index 6ae56d2a..f3fa0eaa 100644 --- a/config/ifort_2021_sse4_mpi.cfg +++ b/config/ifort_2021_sse4_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_xHost.cfg b/config/ifort_2021_xHost.cfg index 1e76a69d..1161833b 100644 --- a/config/ifort_2021_xHost.cfg +++ b/config/ifort_2021_xHost.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=64 -DINTEL diff --git a/configure b/configure index e70820fe..5c38b9f2 100755 --- a/configure +++ b/configure @@ -281,8 +281,8 @@ EOF execute << EOF cd "\${QP_ROOT}"/external - tar -zxf qp2-dependencies/bats-v1.1.0.tar.gz - ( cd bats-core-1.1.0/ ; ./install.sh \${QP_ROOT}) + tar -zxf qp2-dependencies/bats-v1.7.0.tar.gz + ( cd bats-core-1.7.0/ ; ./install.sh \${QP_ROOT}) EOF else diff --git a/etc/qp.rc b/etc/qp.rc index c56661c7..064ca3f7 100644 --- a/etc/qp.rc +++ b/etc/qp.rc @@ -80,8 +80,6 @@ function qp() if [[ -d $NAME ]] ; then [[ -d $EZFIO_FILE ]] && ezfio unset_file ezfio set_file $NAME - else - qp_create_ezfio -h | more fi unset _ARGS ;; diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 90ee61f5..242151e0 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0 +Subproject commit 242151e03d1d6bf042387226431d82d35845686a diff --git a/ocaml/Command_line.ml b/ocaml/Command_line.ml index 1dd57892..602315c6 100644 --- a/ocaml/Command_line.ml +++ b/ocaml/Command_line.ml @@ -1,3 +1,5 @@ +exception Error of string + type short_opt = char type long_opt = string type optional = Mandatory | Optional @@ -181,15 +183,16 @@ let set_specs specs_in = Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]); if show_help () then - (help () ; exit 0); + help () + else + (* Check that all mandatory arguments are set *) + List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs + |> List.iter (fun x -> + match get x.long with + | Some _ -> () + | None -> raise (Error ("--"^x.long^" option is missing.")) + ) - (* Check that all mandatory arguments are set *) - List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs - |> List.iter (fun x -> - match get x.long with - | Some _ -> () - | None -> failwith ("Error: --"^x.long^" option is missing.") - ) ;; diff --git a/ocaml/Command_line.mli b/ocaml/Command_line.mli index 9f6e7022..5ad4ee08 100644 --- a/ocaml/Command_line.mli +++ b/ocaml/Command_line.mli @@ -59,6 +59,8 @@ let () = *) +exception Error of string + type short_opt = char type long_opt = string diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 9b01ac3a..603244c8 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -101,7 +101,7 @@ let to_string_general ~f m = |> String.concat "\n" let to_string = - to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x) + to_string_general ~f:(fun x -> Atom.to_string ~units:Units.Angstrom x) let to_xyz = to_string_general ~f:Atom.to_xyz @@ -113,7 +113,7 @@ let of_xyz_string s = let l = String_ext.split s ~on:'\n' |> List.filter (fun x -> x <> "") - |> list_map (fun x -> Atom.of_string units x) + |> list_map (fun x -> Atom.of_string ~units x) in let ne = ( get_charge { nuclei=l ; diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 270e069f..752a65a0 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -56,3 +56,7 @@ let string_of_string s = s let list_map f l = List.rev_map f l |> List.rev + +let socket_convert socket = + ((Obj.magic (Obj.repr socket)) : [ `Xsub ] Zmq.Socket.t ) + diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index be6c305b..4583b118 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -677,6 +677,7 @@ let run ?o b au c d m p cart xyz_file = let () = + try ( let open Command_line in begin @@ -734,7 +735,7 @@ If a file with the same name as the basis set exists, this file will be read. O let basis = match Command_line.get "basis" with - | None -> assert false + | None -> "" | Some x -> x in @@ -773,10 +774,14 @@ If a file with the same name as the basis set exists, this file will be read. O let xyz_filename = match Command_line.anon_args () with - | [x] -> x - | _ -> (Command_line.help () ; failwith "input file is missing") + | [] -> failwith "input file is missing" + | x::_ -> x in run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename + ) + with + | Failure txt -> Printf.eprintf "Fatal error: %s\n%!" txt + | Command_line.Error txt -> Printf.eprintf "Command line error: %s\n%!" txt diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index d096b15b..dfbab167 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -110,7 +110,7 @@ let run slave ?prefix exe ezfio_file = let task_thread = let thread = Thread.create ( fun () -> - TaskServer.run port_number ) + TaskServer.run ~port:port_number ) in thread (); in diff --git a/ocaml/qp_tunnel.ml b/ocaml/qp_tunnel.ml index 84e50eb5..6885db73 100644 --- a/ocaml/qp_tunnel.ml +++ b/ocaml/qp_tunnel.ml @@ -2,7 +2,7 @@ open Qputils open Qptypes type ezfio_or_address = EZFIO of string | ADDRESS of string -type req_or_sub = REQ | SUB +type req_or_sub = REQ | SUB let localport = 42379 @@ -29,7 +29,7 @@ let () = end; let arg = - let x = + let x = match Command_line.anon_args () with | [x] -> x | _ -> begin @@ -44,7 +44,7 @@ let () = in - let localhost = + let localhost = Lazy.force TaskServer.ip_address in @@ -52,28 +52,28 @@ let () = let long_address = match arg with | ADDRESS x -> x - | EZFIO x -> - let ic = + | EZFIO x -> + let ic = Filename.concat (Qpackage.ezfio_work x) "qp_run_address" |> open_in in - let result = + let result = input_line ic |> String.trim in close_in ic; result in - + let protocol, address, port = match String.split_on_char ':' long_address with | t :: a :: p :: [] -> t, a, int_of_string p - | _ -> failwith @@ + | _ -> failwith @@ Printf.sprintf "%s : Malformed address" long_address in - let zmq_context = + let zmq_context = Zmq.Context.create () in @@ -105,10 +105,10 @@ let () = let create_socket sock_type bind_or_connect addr = - let socket = + let socket = Zmq.Socket.create zmq_context sock_type in - let () = + let () = try bind_or_connect socket addr with @@ -131,37 +131,64 @@ let () = Sys.set_signal Sys.sigint handler; - let new_thread req_or_sub addr_in addr_out = + let new_thread_req addr_in addr_out = let socket_in, socket_out = - match req_or_sub with - | REQ -> create_socket Zmq.Socket.router Zmq.Socket.bind addr_in, create_socket Zmq.Socket.dealer Zmq.Socket.connect addr_out - | SUB -> - create_socket Zmq.Socket.sub Zmq.Socket.connect addr_in, - create_socket Zmq.Socket.pub Zmq.Socket.bind addr_out in - if req_or_sub = SUB then - Zmq.Socket.subscribe socket_in ""; - - - let action_in = - match req_or_sub with - | REQ -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) - | SUB -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + let action_in = + fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out in - let action_out = - match req_or_sub with - | REQ -> (fun () -> Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in ) - | SUB -> (fun () -> () ) + let action_out = + fun () -> Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in in let pollitem = Zmq.Poll.mask_of - [| (socket_in, Zmq.Poll.In) ; (socket_out, Zmq.Poll.In) |] + [| (socket_convert socket_in, Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |] + in + + while !run_status do + + let polling = + Zmq.Poll.poll ~timeout:1000 pollitem + in + + match polling with + | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |] -> ( action_out () ; action_in () ) + | [| _ ; Some Zmq.Poll.In |] -> action_out () + | [| Some Zmq.Poll.In ; _ |] -> action_in () + | _ -> () + done; + + Zmq.Socket.close socket_in; + Zmq.Socket.close socket_out; + in + + let new_thread_sub addr_in addr_out = + let socket_in, socket_out = + create_socket Zmq.Socket.sub Zmq.Socket.connect addr_in, + create_socket Zmq.Socket.pub Zmq.Socket.bind addr_out + in + + Zmq.Socket.subscribe socket_in ""; + + + + let action_in = + fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out + in + + let action_out = + fun () -> () + in + + let pollitem = + Zmq.Poll.mask_of + [| (socket_convert socket_in, Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |] in @@ -173,8 +200,8 @@ let () = match polling with | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |] -> ( action_out () ; action_in () ) - | [| _ ; Some Zmq.Poll.In |] -> action_out () - | [| Some Zmq.Poll.In ; _ |] -> action_in () + | [| _ ; Some Zmq.Poll.In |] -> action_out () + | [| Some Zmq.Poll.In ; _ |] -> action_in () | _ -> () done; @@ -193,8 +220,8 @@ let () = Printf.sprintf "tcp://*:%d" localport in - let f () = - new_thread REQ addr_in addr_out + let f () = + new_thread_req addr_in addr_out in (Thread.create f) () @@ -211,8 +238,8 @@ let () = Printf.sprintf "tcp://*:%d" (localport+2) in - let f () = - new_thread REQ addr_in addr_out + let f () = + new_thread_req addr_in addr_out in (Thread.create f) () in @@ -227,8 +254,8 @@ let () = Printf.sprintf "tcp://*:%d" (localport+1) in - let f () = - new_thread SUB addr_in addr_out + let f () = + new_thread_sub addr_in addr_out in (Thread.create f) () in @@ -236,7 +263,7 @@ let () = let input_thread = - let f () = + let f () = let addr_out = match arg with | EZFIO _ -> None @@ -248,22 +275,22 @@ let () = Printf.sprintf "tcp://*:%d" (localport+9) in - let socket_in = + let socket_in = create_socket Zmq.Socket.rep Zmq.Socket.bind addr_in in let socket_out = - match addr_out with + match addr_out with | Some addr_out -> Some ( create_socket Zmq.Socket.req Zmq.Socket.connect addr_out) | None -> None in - let temp_file = + let temp_file = Filename.temp_file "qp_tunnel" ".tar.gz" in - let get_ezfio_filename () = + let get_ezfio_filename () = match arg with | EZFIO x -> x | ADDRESS _ -> @@ -277,9 +304,9 @@ let () = end in - let get_input () = + let get_input () = match arg with - | EZFIO x -> + | EZFIO x -> begin Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x |> Sys.command |> ignore; @@ -291,11 +318,11 @@ let () = in ignore @@ Unix.lseek fd 0 Unix.SEEK_SET ; let bstr = - Unix.map_file fd Bigarray.char + Unix.map_file fd Bigarray.char Bigarray.c_layout false [| len |] |> Bigarray.array1_of_genarray in - let result = + let result = String.init len (fun i -> bstr.{i}) ; in Unix.close fd; @@ -313,7 +340,7 @@ let () = end in - let () = + let () = match socket_out with | None -> () | Some socket_out -> @@ -329,7 +356,7 @@ let () = | ADDRESS _ -> begin Printf.printf "Getting input... %!"; - let ezfio_filename = + let ezfio_filename = get_ezfio_filename () in Printf.printf "%s%!" ezfio_filename; @@ -343,7 +370,7 @@ let () = |> Sys.command |> ignore ; let oc = Filename.concat (Qpackage.ezfio_work ezfio_filename) "qp_run_address" - |> open_out + |> open_out in Printf.fprintf oc "tcp://%s:%d\n" localhost localport; close_out oc; @@ -359,9 +386,9 @@ let () = let action () = match Zmq.Socket.recv socket_in with | "get_input" -> get_input () - |> Zmq.Socket.send socket_in + |> Zmq.Socket.send socket_in | "get_ezfio_filename" -> get_ezfio_filename () - |> Zmq.Socket.send socket_in + |> Zmq.Socket.send socket_in | "test" -> Zmq.Socket.send socket_in "OK" | x -> Printf.sprintf "Message '%s' not understood" x |> Zmq.Socket.send socket_in @@ -372,7 +399,7 @@ On remote hosts, create ssh tunnel using: ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s & Or from this host connect to clients using: ssh -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d & -%!" +%!" (port ) localhost (localport ) (port+1) localhost (localport+1) (port+2) localhost (localport+2) @@ -392,12 +419,12 @@ Or from this host connect to clients using: match polling.(0) with | Some Zmq.Poll.In -> action () | None -> () - | Some Zmq.Poll.In_out + | Some Zmq.Poll.In_out | Some Zmq.Poll.Out -> () done; - let () = + let () = match socket_out with | Some socket_out -> Zmq.Socket.close socket_out | None -> () @@ -415,7 +442,7 @@ Or from this host connect to clients using: Thread.join ocaml_thread; Zmq.Context.terminate zmq_context; Printf.printf "qp_tunnel exited properly.\n" - + diff --git a/scripts/cipsi_save.sh b/scripts/cipsi_save.sh new file mode 100644 index 00000000..a4d9b65e --- /dev/null +++ b/scripts/cipsi_save.sh @@ -0,0 +1,27 @@ +#!/bin/bash +# +# This script runs a CIPSI calculation as a sequence of single CIPSI iterations. +# After each iteration, the EZFIO directory is saved. +# +# Usage: cipsi_save [EZFIO_FILE] [NDET] +# +# Example: cipsi_save file.ezfio 10000 + +EZ=$1 +NDETMAX=$2 + +qp set_file ${EZ} +qp reset -d +qp set determinants read_wf true +declare -i NDET +NDET=1 +while [[ ${NDET} -lt ${NDETMAX} ]] +do + NDET=$(($NDET + $NDET)) + qp set determinants n_det_max $NDET + qp run fci > ${EZ}.out + NDET=$(qp get determinants n_det) + mv ${EZ}.out ${EZ}.${NDET}.out + cp -r ${EZ} ${EZ}.${NDET} +done + diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 7fcf19eb..e01359c5 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -34,3 +34,9 @@ doc: Maximum number of excitation for beta determinants with respect to the Hart interface: ezfio,ocaml,provider default: -1 +[twice_hierarchy_max] +type: integer +doc: Twice the maximum hierarchy parameter (excitation degree plus half the seniority number). Using -1 selects all determinants +interface: ezfio,ocaml,provider +default: -1 + diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index c7cee1ac..1328e7a0 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -290,9 +290,13 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call set_multiple_levels_omp(.False.) - print '(A)', '========== ======================= ===================== ===================== ===========' - print '(A)', ' Samples Energy Variance Norm^2 Seconds' - print '(A)', '========== ======================= ===================== ===================== ===========' + ! old + !print '(A)', '========== ======================= ===================== ===================== ===========' + !print '(A)', ' Samples Energy Variance Norm^2 Seconds' + !print '(A)', '========== ======================= ===================== ===================== ===========' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + print '(A)', ' Samples Energy PT2 Variance Norm^2 Convergence Seconds' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' PROVIDE global_selection_buffer @@ -316,7 +320,10 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') call set_multiple_levels_omp(.True.) - print '(A)', '========== ======================= ===================== ===================== ===========' + ! old + !print '(A)', '========== ======================= ===================== ===================== ===========' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) @@ -414,6 +421,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ double precision :: rss double precision, external :: memory_of_double, memory_of_int + character(len=20) :: format_str1, str_error1, format_str2, str_error2 + character(len=20) :: format_str3, str_error3, format_str4, str_error4 + character(len=20) :: format_value1, format_value2, format_value3, format_value4 + character(len=20) :: str_value1, str_value2, str_value3, str_value4 + character(len=20) :: str_conv + double precision :: value1, value2, value3, value4 + double precision :: error1, error2, error3, error4 + integer :: size1,size2,size3,size4 + + double precision :: conv_crit + sending =.False. rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2) @@ -537,14 +555,60 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, & - pt2_data % pt2(pt2_stoch_istate) +E, & - pt2_data_err % pt2(pt2_stoch_istate), & - pt2_data % variance(pt2_stoch_istate), & - pt2_data_err % variance(pt2_stoch_istate), & - pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & - pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & - time-time0 + + value1 = pt2_data % pt2(pt2_stoch_istate) + E + error1 = pt2_data_err % pt2(pt2_stoch_istate) + value2 = pt2_data % pt2(pt2_stoch_istate) + error2 = pt2_data_err % pt2(pt2_stoch_istate) + value3 = pt2_data % variance(pt2_stoch_istate) + error3 = pt2_data_err % variance(pt2_stoch_istate) + value4 = pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate) + error4 = pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate) + + ! Max size of the values (FX.Y) with X=size + size1 = 15 + size2 = 12 + size3 = 12 + size4 = 12 + + ! To generate the format: number(error) + call format_w_error(value1,error1,size1,8,format_value1,str_error1) + call format_w_error(value2,error2,size2,8,format_value2,str_error2) + call format_w_error(value3,error3,size3,8,format_value3,str_error3) + call format_w_error(value4,error4,size4,8,format_value4,str_error4) + + ! value > string with the right format + write(str_value1,'('//format_value1//')') value1 + write(str_value2,'('//format_value2//')') value2 + write(str_value3,'('//format_value3//')') value3 + write(str_value4,'('//format_value4//')') value4 + + ! Convergence criterion + conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & + (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) + write(str_conv,'(G10.3)') conv_crit + + write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,& + adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),& + adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),& + adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),& + adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),& + adjustl(str_conv),& + time-time0 + + ! Old print + !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, & + ! pt2_data % pt2(pt2_stoch_istate) +E, & + ! pt2_data_err % pt2(pt2_stoch_istate), & + ! pt2_data % variance(pt2_stoch_istate), & + ! pt2_data_err % variance(pt2_stoch_istate), & + ! pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & + ! pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & + ! time-time0, & + ! pt2_data % pt2(pt2_stoch_istate), & + ! dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & + ! (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) + if (stop_now .or. ( & (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 91bd3a38..de7c209c 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -61,10 +61,14 @@ subroutine run_selection_slave(thread,iproc,energy) if (N /= buf%N) then print *, 'N=', N print *, 'buf%N=', buf%N - print *, 'bug in ', irp_here - stop '-1' + print *, 'In ', irp_here, ': N /= buf%N' + stop -1 end if end if + if (i_generator > N_det_generators) then + print *, 'In ', irp_here, ': i_generator > N_det_generators' + stop -1 + endif call select_connected(i_generator,energy,pt2_data,buf,subset,pt2_F(i_generator)) endif diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index d4d44c2d..ec60c606 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -195,7 +195,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer :: l_a, nmax, idx integer, allocatable :: indices(:), exc_degree(:), iorder(:) - double precision, parameter :: norm_thr = 1.d-16 + + ! Removed to avoid introducing determinants already presents in the wf + !double precision, parameter :: norm_thr = 1.d-16 + allocate (indices(N_det), & exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) @@ -215,10 +218,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d i = psi_bilinear_matrix_rows(l_a) if (nt + exc_degree(i) <= 4) then idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -242,10 +246,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d idx = psi_det_sorted_order( & psi_bilinear_matrix_order( & psi_bilinear_matrix_transp_order(l_a))) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -253,8 +258,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d deallocate(exc_degree) nmax=k-1 - call isort_noidx(indices,nmax) - ! Start with 32 elements. Size will double along with the filtering. allocate(preinteresting(0:32), prefullinteresting(0:32), & interesting(0:32), fullinteresting(0:32)) @@ -566,6 +569,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision, external :: diag_H_mat_elem_fock double precision :: E_shift double precision :: s_weight(N_states,N_states) + logical, external :: is_in_wavefunction PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs do jstate=1,N_states do istate=1,N_states @@ -707,6 +711,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (do_cycle) cycle endif + if (twice_hierarchy_max >= 0) then + s = 0 + do k=1,N_int + s = s + popcnt(ieor(det(k,1),det(k,2))) + enddo + if ( mod(s,2)>0 ) stop 'For now, hierarchy CI is defined only for an even number of electrons' + if (excitation_ref == 1) then + call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int) + else if (excitation_ref == 2) then + stop 'For now, hierarchy CI is defined only for a single reference determinant' +! do k=1,N_dominant_dets_of_cfgs +! call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) +! enddo + endif + integer :: twice_hierarchy + twice_hierarchy = degree + s/2 + if (twice_hierarchy > twice_hierarchy_max) cycle + endif + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) w = 0d0 @@ -777,7 +800,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi = mat(istate, p1, p2) - pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate) + do k=1,N_states + pt2_data % overlap(k,istate) = pt2_data % overlap(k,istate) + coef(k) * coef(istate) + end do pt2_data % variance(istate) = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate) @@ -828,8 +853,27 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif end select + + ! To force the inclusion of determinants with a positive pt2 contribution + if (e_pert(istate) > 1d-8) then + w = -huge(1.0) + endif + end do +!!!BEGIN_DEBUG +! ! To check if the pt2 is taking determinants already in the wf +! if (is_in_wavefunction(det(N_int,1),N_int)) then +! print*, 'A determinant contributing to the pt2 is already in' +! print*, 'the wave function:' +! call print_det(det(N_int,1),N_int) +! print*,'contribution to the pt2 for the states:', e_pert(:) +! print*,'error in the filtering in' +! print*, 'cipsi/selection.irp.f sub: selecte_singles_and_doubles' +! print*, 'abort' +! call abort +! endif +!!!END_DEBUG integer(bit_kind) :: occ(N_int,2), n if (h0_type == 'CFG') then diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index f72197c2..2b16a5f7 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -74,7 +74,7 @@ subroutine run print*,'******************************************************' print*,'Excitation energies (au) (eV)' do i = 2, N_states - print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0 + print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev enddo print*,'' endif diff --git a/src/cis_read/EZFIO.cfg b/src/cis_read/EZFIO.cfg new file mode 100644 index 00000000..955d1bef --- /dev/null +++ b/src/cis_read/EZFIO.cfg @@ -0,0 +1,7 @@ +[energy] +type: double precision +doc: Variational |CIS| energy +interface: ezfio +size: (determinants.n_states) + + diff --git a/src/cis_read/NEED b/src/cis_read/NEED new file mode 100644 index 00000000..42992ac6 --- /dev/null +++ b/src/cis_read/NEED @@ -0,0 +1,3 @@ +selectors_full +generators_full +davidson_undressed diff --git a/src/cis_read/README.rst b/src/cis_read/README.rst new file mode 100644 index 00000000..31648636 --- /dev/null +++ b/src/cis_read/README.rst @@ -0,0 +1,5 @@ +=== +cis_read +=== + +Reads the input WF and performs all singles on top of it. diff --git a/src/cis_read/cis_read.irp.f b/src/cis_read/cis_read.irp.f new file mode 100644 index 00000000..055b5e15 --- /dev/null +++ b/src/cis_read/cis_read.irp.f @@ -0,0 +1,88 @@ +program cis + implicit none + BEGIN_DOC +! +! Configuration Interaction with Single excitations. +! +! This program takes a reference Slater determinant of ROHF-like +! occupancy, and performs all single excitations on top of it. +! Disregarding spatial symmetry, it computes the `n_states` lowest +! eigenstates of that CI matrix. (see :option:`determinants n_states`) +! +! This program can be useful in many cases: +! +! +! 1. Ground state calculation +! +! To be sure to have the lowest |SCF| solution, perform an :ref:`scf` +! (see the :ref:`module_hartree_fock` module), then a :ref:`cis`, save the +! natural orbitals (see :ref:`save_natorb`) and re-run an :ref:`scf` +! optimization from this |MO| guess. +! +! +! 2. Excited states calculations +! +! The lowest excited states are much likely to be dominated by +! single-excitations. Therefore, running a :ref:`cis` will save the +! `n_states` lowest states within the |CIS| space in the |EZFIO| +! directory, which can afterwards be used as guess wave functions for +! a further multi-state |FCI| calculation if :option:`determinants +! read_wf` is set to |true| before running the :ref:`fci` executable. +! +! +! If :option:`determinants s2_eig` is set to |true|, the |CIS| +! will only retain states having the expected |S^2| value (see +! :option:`determinants expected_s2`). Otherwise, the |CIS| will take +! the lowest :option:`determinants n_states`, whatever multiplicity +! they are. +! +! .. note:: +! +! To discard some orbitals, use the :ref:`qp_set_mo_class` +! command to specify: +! +! * *core* orbitals which will be always doubly occupied +! +! * *act* orbitals where an electron can be either excited from or to +! +! * *del* orbitals which will be never occupied +! + END_DOC + read_wf = .True. + TOUCH read_wf + call run +end + +subroutine run + implicit none + integer :: i + + + if(pseudo_sym)then + call H_apply_cis_sym + else + call H_apply_cis + endif + print*,'' + 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*,'******************************************************' + print*,'Excitation energies (au) (eV)' + do i = 2, N_states + print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0 + enddo + print*,'' + endif + + call ezfio_set_cis_energy(CI_energy) + psi_coef = ci_eigenvectors + SOFT_TOUCH psi_coef + call save_wavefunction_truncated(save_threshold) + +end diff --git a/src/cis_read/h_apply.irp.f b/src/cis_read/h_apply.irp.f new file mode 100644 index 00000000..14389bed --- /dev/null +++ b/src/cis_read/h_apply.irp.f @@ -0,0 +1,14 @@ +! Generates subroutine H_apply_cis +! -------------------------------- + +BEGIN_SHELL [ /usr/bin/env python3 ] +from generate_h_apply import H_apply +H = H_apply("cis",do_double_exc=False) +print(H) + +H = H_apply("cis_sym",do_double_exc=False) +H.filter_only_connected_to_hf() +print(H) + +END_SHELL + diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats index 58d996f8..6e110aa3 100644 --- a/src/cisd/30.cisd.bats +++ b/src/cisd/30.cisd.bats @@ -77,7 +77,7 @@ function run() { [[ -n $TRAVIS ]] && skip qp set_file ch4.ezfio qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" - run -40.2403962667047 -39.8433221754964 + run -40.2403962667047 -39.843315 } @test "SiH3" { # 20.2202s 1.38648m diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f index fca3b10e..3e1e8d97 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -47,6 +47,37 @@ program cisd PROVIDE N_states read_wf = .False. SOFT_TOUCH read_wf + + integer :: i,k + + if(pseudo_sym)then + call H_apply_cisd_sym + else + call H_apply_cisd + endif + double precision :: r1, r2 + double precision, allocatable :: U_csf(:,:) + + allocate(U_csf(N_csf,N_states)) + U_csf = 0.d0 + U_csf(1,1) = 1.d0 + do k=2,N_states + do i=1,N_csf + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dacos(-1.d0)*2.d0*r2 + U_csf(i,k) = r1*dcos(r2) + enddo + U_csf(k,k) = U_csf(k,k) +100.d0 + enddo + do k=1,N_states + call normalize(U_csf(1,k),N_csf) + enddo + call convertWFfromCSFtoDET(N_states,U_csf(1,1),psi_coef(1,1)) + deallocate(U_csf) + SOFT_TOUCH psi_coef + call run end @@ -56,20 +87,16 @@ subroutine run 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_truncated(save_threshold) 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) + if (elec_alpha_num + elec_beta_num >= 4) then + cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2) + endif enddo print *, 'N_det = ', N_det print*,'' @@ -78,26 +105,43 @@ subroutine run 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 (elec_alpha_num + elec_beta_num >= 4) then + print*,'' + print*,'******************************' + print *, 'CISD+Q Energies' + do i = 1,N_states + print *, i, cisdq(i) + enddo + endif 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 + if (elec_alpha_num + elec_beta_num >= 4) 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)) * ha_to_ev, & + (cisdq(i) - cisdq(1)) * ha_to_ev + enddo + else + print*,'' + print*,'******************************' + print*,'Excitation energies (au) (CISD)' + do i = 2, N_states + print*, i ,CI_energy(i) - CI_energy(1) + enddo + print*,'' + print*,'******************************' + print*,'Excitation energies (eV) (CISD)' + do i = 2, N_states + print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev + enddo + endif endif end diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index 23b984a0..746de04e 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -68,10 +68,16 @@ void getBFIndexList(int NSOMO, int *BF1, int *IdxListBF1){ break; } } - BFcopy[Iidx] = -1; - BFcopy[Jidx] = -1; - IdxListBF1[Jidx] = Iidx; - IdxListBF1[Iidx] = Jidx; + if(countN1 <= 0){ + BFcopy[Iidx] = -1; + IdxListBF1[Iidx] = Iidx; + } + else{ + BFcopy[Iidx] = -1; + BFcopy[Jidx] = -1; + IdxListBF1[Jidx] = Iidx; + IdxListBF1[Iidx] = Jidx; + } } } @@ -328,10 +334,21 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl Get Overlap ************************************/ // Fill matrix + + int rowsbftodetI, colsbftodetI; + + /*********************************** + Get BFtoDeterminant Matrix + ************************************/ + + printf("In convertcsftodet\n"); + convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); + int rowsI = 0; int colsI = 0; - getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); + //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); + getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); /*********************************** @@ -342,14 +359,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI); - /*********************************** - Get BFtoDeterminant Matrix - ************************************/ - - int rowsbftodetI, colsbftodetI; - - convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); - /*********************************** Get Final CSF to Det Matrix ************************************/ @@ -1297,12 +1306,16 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv double fac = 1.0; for(int i = 0; i < NSOMO; i++) donepq[i] = 0.0; + for(int i=0;i 0.0 || donepq[idxq] > 0.0) continue; + if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0 || idxp == idxq) continue; fac *= 2.0; donepq[idxp] = 1.0; donepq[idxq] = 1.0; @@ -1319,15 +1332,19 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv } shft /= 2; } - + // Now get the addresses int inpdet[NSOMO]; int phase_cfg_to_qp=1; int addr = -1; for(int i = 0; i < npairs; i++){ - for(int j = 0; j < NSOMO; j++) + for(int j = 0; j < NSOMO; j++) { inpdet[j] = detslist[i*NSOMO + j]; + printf(" %d ",inpdet[j]); + } + printf("\n"); findAddofDetDriver(dettree, NSOMO, inpdet, &addr); + printf("(%d) - addr = %d\n",i,addr); // Calculate the phase for cfg to QP2 conversion //get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp); //rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac); @@ -1416,7 +1433,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int * getIthBFDriver(&bftree, NSOMO, addI, BF1); getBFIndexList(NSOMO, BF1, IdxListBF1); - // Get ith row getbftodetfunction(&dettree, NSOMO, MS, IdxListBF1, rowvec); @@ -1425,6 +1441,11 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int * for(int k=0;k NSOMOMax || icpl < 0 || izeros > zeromax ) return; // If we find a valid BF assign its address - if(isomo == NSOMOMax){ + if(isomo == NSOMOMax && icpl == MSmax){ (*inode)->addr = bftree->rootNode->addr; bftree->rootNode->addr += 1; return; } // Call 0 branch - if(((*inode)->C0) == NULL && izeros+1 <= zeromax){ - ((*inode)->C0) = malloc(sizeof(Node)); - (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo }; - buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax); + if(izeros+1 <= zeromax){ + if(((*inode)->C0) == NULL){ + ((*inode)->C0) = malloc(sizeof(Node)); + (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo }; + buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax); + } + else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax); } - else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax); // Call 1 branch - if(((*inode)->C1) == NULL && icpl-1 >= 0){ - ((*inode)->C1) = malloc(sizeof(Node)); - (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo }; - buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax); + if(icpl-1 >=0){ + if(((*inode)->C1) == NULL){ + ((*inode)->C1) = malloc(sizeof(Node)); + (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo }; + buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax); + } + else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax); } - else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax); return; } diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 0c3c6f92..7aaaa842 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -124,7 +124,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N stop -1 endif - itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itermax = max(2,min(davidson_sze_max, sze_csf/N_st_diag))+1 itertot = 0 if (state_following) then @@ -263,29 +263,20 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! =================== converged = .False. - + call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),U_csf(1,1)) do k=N_st+1,N_st_diag - do i=1,sze + do i=1,sze_csf call random_number(r1) call random_number(r2) r1 = dsqrt(-2.d0*dlog(r1)) r2 = dtwo_pi*r2 - u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + U_csf(i,k) = r1*dcos(r2) * u_csf(i,k-N_st) enddo - u_in(k,k) = u_in(k,k) + 10.d0 + U_csf(k,k) = u_csf(k,k) + 10.d0 enddo do k=1,N_st_diag - call normalize(u_in(1,k),sze) + call normalize(U_csf(1,k),sze_csf) enddo - - do k=1,N_st_diag - do i=1,sze - U(i,k) = u_in(i,k) - enddo - enddo - - ! Make random verctors eigenstates of S2 - call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1)) call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1)) do while (.not.converged) diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 46ad8f78..6930cc07 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -1,3 +1,13 @@ +BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ] + implicit none + BEGIN_DOC + ! If 'det', use in Davidson + ! + ! If 'cfg', use in Davidson + END_DOC + sigma_vector_algorithm = 'det' +END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] implicit none BEGIN_DOC @@ -61,9 +71,18 @@ END_PROVIDER if (diag_algorithm == "Davidson") then if (do_csf) then - call davidson_diag_H_csf(psi_det,CI_eigenvectors, & - size(CI_eigenvectors,1),CI_electronic_energy, & - N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged) + if (sigma_vector_algorithm == 'det') then + call davidson_diag_H_csf(psi_det,CI_eigenvectors, & + size(CI_eigenvectors,1),CI_electronic_energy, & + N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged) +! else if (sigma_vector_algorithm == 'cfg') then +! call davidson_diag_H_csf(psi_det,CI_eigenvectors, & +! size(CI_eigenvectors,1),CI_electronic_energy, & +! N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged) +! else +! print *, irp_here +! stop 'bug' + endif else call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_s2, & size(CI_eigenvectors,1),CI_electronic_energy, & diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index 5b12a6d9..52092fad 100644 --- a/src/determinants/determinants.irp.f +++ b/src/determinants/determinants.irp.f @@ -77,28 +77,31 @@ BEGIN_PROVIDER [ integer, psi_det_size ] END_DOC PROVIDE ezfio_filename logical :: exists - if (mpi_master) then - call ezfio_has_determinants_n_det(exists) - if (exists) then - call ezfio_get_determinants_n_det(psi_det_size) - else - psi_det_size = 1 + psi_det_size = N_states + PROVIDE mpi_master + if (read_wf) then + if (mpi_master) then + call ezfio_has_determinants_n_det(exists) + if (exists) then + call ezfio_get_determinants_n_det(psi_det_size) + else + psi_det_size = N_states + endif + call write_int(6,psi_det_size,'Dimension of the psi arrays') endif - call write_int(6,psi_det_size,'Dimension of the psi arrays') + 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( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read psi_det_size with MPI' + endif + IRP_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( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read psi_det_size with MPI' - endif - IRP_ENDIF - END_PROVIDER @@ -539,7 +542,7 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) integer :: i,j,k, ndet_qp_edit if (mpi_master) then - ndet_qp_edit = min(ndet,N_det_qp_edit) + ndet_qp_edit = min(ndet,10000) call ezfio_set_determinants_N_int(N_int) call ezfio_set_determinants_bit_kind(bit_kind) diff --git a/src/determinants/dipole_moments.irp.f b/src/determinants/dipole_moments.irp.f index 06fca0cd..b411dda4 100644 --- a/src/determinants/dipole_moments.irp.f +++ b/src/determinants/dipole_moments.irp.f @@ -66,9 +66,13 @@ END_PROVIDER write(*,'(i16)',advance='no') i end do write(*,*) '' - write(*,'(A17,100(1pE16.8))') 'x_dipole_moment = ',x_dipole_moment - write(*,'(A17,100(1pE16.8))') 'y_dipole_moment = ',y_dipole_moment - write(*,'(A17,100(1pE16.8))') 'z_dipole_moment = ',z_dipole_moment + write(*,'(A23,100(1pE16.8))') 'x_dipole_moment (au) = ',x_dipole_moment + write(*,'(A23,100(1pE16.8))') 'y_dipole_moment (au) = ',y_dipole_moment + write(*,'(A23,100(1pE16.8))') 'z_dipole_moment (au) = ',z_dipole_moment + write(*,*) '' + write(*,'(A23,100(1pE16.8))') 'x_dipole_moment (D) = ',x_dipole_moment * au_to_D + write(*,'(A23,100(1pE16.8))') 'y_dipole_moment (D) = ',y_dipole_moment * au_to_D + write(*,'(A23,100(1pE16.8))') 'z_dipole_moment (D) = ',z_dipole_moment * au_to_D !print*, 'x_dipole_moment = ',x_dipole_moment !print*, 'y_dipole_moment = ',y_dipole_moment !print*, 'z_dipole_moment = ',z_dipole_moment diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 3a33a37d..897607a9 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -623,7 +623,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase integer :: n_occ_ab(2) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE ao_one_e_integrals mo_one_e_integrals ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -681,7 +682,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij) case (1) call get_single_excitation(key_i,key_j,exc,phase,Nint) !DIR$ FORCEINLINE - call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) if (exc(0,1,1) == 1) then ! Single alpha m = exc(1,1,1) @@ -700,10 +700,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij) end select end - - - - subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase) use bitmasks implicit none @@ -1038,7 +1034,6 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) end - subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none diff --git a/src/determinants/slater_rules_wee_mono.irp.f b/src/determinants/slater_rules_wee_mono.irp.f index 4c1c9330..7c2ad148 100644 --- a/src/determinants/slater_rules_wee_mono.irp.f +++ b/src/determinants/slater_rules_wee_mono.irp.f @@ -282,9 +282,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij) double precision :: get_two_e_integral integer :: m,n,p,q integer :: i,j,k - integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase,phase_2 - integer :: n_occ_ab(2) PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy ASSERT (Nint > 0) @@ -342,7 +340,6 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij) case (1) call get_single_excitation(key_i,key_j,exc,phase,Nint) !DIR$ FORCEINLINE - call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) if (exc(0,1,1) == 1) then ! Mono alpha m = exc(1,1,1) diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index f2c78067..ac34251f 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -8,12 +8,12 @@ function run() { test_exe fci || skip qp edit --check qp set perturbation do_pt2 False - qp set determinants n_det_max 8000 + 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 8 qp run fci - energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)" + energy1="$(qp get fci energy | tr '[]' ' ' | cut -d ',' -f 1)" eq $energy1 $1 $thresh } @@ -22,166 +22,167 @@ function run_stoch() { thresh=$2 test_exe fci || skip qp set perturbation do_pt2 True + qp set perturbation pt2_relative_error 0.005 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 qp run fci - energy1="$(ezfio get fci energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)" + energy1="$(qp get fci energy_extrapolated | tr '[]' ' ' | cut -d ',' -f 1)" eq $energy1 $1 $thresh } -@test "B-B" { + +@test "B-B" { # 0:00:10 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 + run_stoch -49.14097596 0.0001 10000 } -@test "F2" { # 4.07m - [[ -n $TRAVIS ]] && skip - qp set_file f2.ezfio - qp set_frozen_core - run_stoch -199.304922384814 3.e-4 100000 -} - -@test "NH3" { # 10.6657s +@test "NH3" { # 0:00:11 qp set_file nh3.ezfio qp set_mo_class --core="[1-4]" --act="[5-72]" - run -56.244753429144986 3.e-4 100000 + run -56.24474908 1.e-5 10000 } -@test "DHNO" { # 11.4721s +@test "DHNO" { # 0:00:10 qp set_file dhno.ezfio qp set_mo_class --core="[1-7]" --act="[8-64]" - run -130.459020029816 3.e-4 100000 + run -130.45904647 1.e-4 10000 } -@test "HCO" { # 12.2868s +@test "HCO" { # 0:01:16 qp set_file hco.ezfio - run -113.389297812482 6.e-4 100000 + run_stoch -113.41448940 2.e-3 50000 } -@test "H2O2" { # 12.9214s +@test "H2O2" { # 0:01:48 qp set_file h2o2.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]" - run -151.00467 1.e-4 100000 + run_stoch -151.02437936 2.e-3 100000 } -@test "HBO" { # 13.3144s +@test "HBO" { # 0:00:46 [[ -n $TRAVIS ]] && skip qp set_file hbo.ezfio - run -100.212560384678 1.e-3 100000 + run_stoch -100.221198108988 2.e-3 50000 } -@test "H2O" { # 11.3727s +@test "H2O" { # 0:01:05 [[ -n $TRAVIS ]] && skip qp set_file h2o.ezfio - run -76.2361605151999 3.e-4 100000 + run_stoch -76.241332121813 1.e-3 100000 } -@test "ClO" { # 13.3755s +@test "ClO" { # 0:03:07 [[ -n $TRAVIS ]] && skip qp set_file clo.ezfio - run -534.545616787223 3.e-4 100000 + run_stoch -534.573564655419 1.e-3 100000 } -@test "SO" { # 13.4952s +@test "SO" { # 0:01:49 [[ -n $TRAVIS ]] && skip qp set_file so.ezfio - run -26.0096209515081 1.e-3 100000 + run_stoch -26.04335528 5.e-3 100000 } -@test "H2S" { # 13.6745s +@test "H2S" { # 0:01:12 [[ -n $TRAVIS ]] && skip qp set_file h2s.ezfio - run -398.859168655255 3.e-4 100000 + run_stoch -398.865173546866 1.e-3 50000 } -@test "OH" { # 13.865s +@test "OH" { # 0:00:41 [[ -n $TRAVIS ]] && skip qp set_file oh.ezfio - run -75.6121856748294 3.e-4 100000 + run_stoch -75.6193013819546 1.e-3 50000 } -@test "SiH2_3B1" { # 13.938ss +@test "SiH2_3B1" { # 0:00:50 [[ -n $TRAVIS ]] && skip qp set_file sih2_3b1.ezfio - run -290.0175411299477 3.e-4 100000 + run_stoch -290.01754869 3.e-5 50000 } -@test "H3COH" { # 14.7299s +@test "H3COH" { # 0:01:05 [[ -n $TRAVIS ]] && skip qp set_file h3coh.ezfio - run -115.205191406072 3.e-4 100000 + run_stoch -115.224147057725 2.e-3 50000 } -@test "SiH3" { # 15.99s +@test "SiH3" { # 0:01:09 [[ -n $TRAVIS ]] && skip qp set_file sih3.ezfio - run -5.57241217753818 3.e-4 100000 + run_stoch -5.57812512359276 1.e-3 50000 } -@test "CH4" { # 16.1612s +@test "CH4" { # 0:02:06 [[ -n $TRAVIS ]] && skip qp set_file ch4.ezfio qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" - run -40.2409678239136 3.e-4 100000 + run_stoch -40.2419474611994 1.e-4 100000 } -@test "ClF" { # 16.8864s +@test "ClF" { # 0:01:55 [[ -n $TRAVIS ]] && skip qp set_file clf.ezfio - run -559.169313755572 3.e-4 100000 + run_stoch -559.20666465 1.e-2 50000 } -@test "SO2" { # 17.5645s +@test "SO2" { # 0:00:24 [[ -n $TRAVIS ]] && skip qp set_file so2.ezfio qp set_mo_class --core="[1-8]" --act="[9-87]" - run -41.5746738713298 3.e-4 100000 + run_stoch -41.57468756 1.e-4 50000 } -@test "C2H2" { # 17.6827s +@test "C2H2" { # 0:00:57 [[ -n $TRAVIS ]] && skip qp set_file c2h2.ezfio qp set_mo_class --act="[1-30]" --del="[31-36]" - run -12.3685464085969 3.e-4 100000 + run_stoch -12.3862664765532 1.e-3 50000 } -@test "N2" { # 18.0198s +@test "N2" { # 0:01:15 [[ -n $TRAVIS ]] && skip qp set_file n2.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]" - run -109.28681540699360 3.e-4 100000 + run_stoch -109.311954243348 2.e-3 50000 } -@test "N2H4" { # 18.5006s +@test "N2H4" { # 0:00:51 [[ -n $TRAVIS ]] && skip qp set_file n2h4.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]" - run -111.367332681559 3.e-4 100000 + run_stoch -111.38119165053 1.e-3 50000 } -@test "CO2" { # 21.1748s +@test "CO2" { # 0:01:00 [[ -n $TRAVIS ]] && skip qp set_file co2.ezfio qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]" - run -187.968547952413 3.e-4 100000 + run_stoch -188.002190327443 2.e-3 50000 } - -@test "[Cu(NH3)4]2+" { # 25.0417s +@test "[Cu(NH3)4]2+" { # 0:01:53 [[ -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.9869374387192 3.e-04 100000 + run_stoch -1862.98705340328 1.e-05 50000 } -@test "HCN" { # 20.3273s +@test "HCN" { # 0:01:26 [[ -n $TRAVIS ]] && skip qp set_file hcn.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]" - run -93.0771143355433 3.e-4 100000 + run_stoch -93.0980746734051 5.e-4 50000 +} + +@test "F2" { # 0:03:34 + [[ -n $TRAVIS ]] && skip + qp set_file f2.ezfio + qp set_frozen_core + run_stoch -199.307512211742 0.002 100000 } diff --git a/src/fci/EZFIO.cfg b/src/fci/EZFIO.cfg index d897428a..6b46a85f 100644 --- a/src/fci/EZFIO.cfg +++ b/src/fci/EZFIO.cfg @@ -10,3 +10,9 @@ doc: Calculated |FCI| energy + |PT2| interface: ezfio size: (determinants.n_states) +[energy_extrapolated] +type: double precision +doc: Calculated |FCI| energy + |PT2| +interface: ezfio +size: (determinants.n_states) + diff --git a/src/iterations/print_extrapolation.irp.f b/src/iterations/print_extrapolation.irp.f index cb46fb67..7c6dbb9b 100644 --- a/src/iterations/print_extrapolation.irp.f +++ b/src/iterations/print_extrapolation.irp.f @@ -35,12 +35,13 @@ subroutine print_extrapolated_energy do k=2,min(N_iter,8) write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), & extrapolated_energy(k,i) - extrapolated_energy(k,1), & - (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0 + (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * ha_to_ev enddo write(*,*) '=========== ', '=================== ', '=================== ', '===================' enddo print *, '' + call ezfio_set_fci_energy_extrapolated(extrapolated_energy(min(N_iter,3),1:N_states)) end subroutine diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 8e6285e2..a0db3534 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -36,7 +36,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s write(*,fmt) '# E ', e_(1:N_states_p) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) - write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 + write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*ha_to_ev endif write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) @@ -47,8 +47,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) - write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & - dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*ha_to_ev, & + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*ha_to_ev, k=1,N_states_p) endif write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) @@ -82,19 +82,19 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s print *, 'Variational Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i) - e_(1)), & - (e_(i) - e_(1)) * 27.211396641308d0 + (e_(i) - e_(1)) * ha_to_ev enddo print *, '-----' print*, 'Variational + perturbative Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), & - (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0 + (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * ha_to_ev enddo print *, '-----' print*, 'Variational + renormalized perturbative Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), & - (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0 + (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * ha_to_ev enddo endif diff --git a/src/mo_two_e_erf_ints/map_integrals_erf.irp.f b/src/mo_two_e_erf_ints/map_integrals_erf.irp.f index 73050ec5..3405ec2b 100644 --- a/src/mo_two_e_erf_ints/map_integrals_erf.irp.f +++ b/src/mo_two_e_erf_ints/map_integrals_erf.irp.f @@ -235,11 +235,11 @@ subroutine get_mo_two_e_integrals_erf_ij(k,l,sze,out_array,map) logical :: integral_is_in_map if (key_kind == 8) then - call i8radix_sort(hash,iorder,kk,-1) + call i8sort(hash,iorder,kk) else if (key_kind == 4) then - call iradix_sort(hash,iorder,kk,-1) + call isort(hash,iorder,kk) else if (key_kind == 2) then - call i2radix_sort(hash,iorder,kk,-1) + call i2sort(hash,iorder,kk) endif call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk) @@ -290,11 +290,11 @@ subroutine get_mo_two_e_integrals_erf_i1j1(k,l,sze,out_array,map) logical :: integral_is_in_map if (key_kind == 8) then - call i8radix_sort(hash,iorder,kk,-1) + call i8sort(hash,iorder,kk) else if (key_kind == 4) then - call iradix_sort(hash,iorder,kk,-1) + call isort(hash,iorder,kk) else if (key_kind == 2) then - call i2radix_sort(hash,iorder,kk,-1) + call i2sort(hash,iorder,kk) endif call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk) diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index d58932ce..6f4c5c17 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -53,7 +53,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] ! call four_idx_novvvv call four_idx_novvvv_old else - call add_integrals_to_map(full_ijkl_bitmask_4) + if (32.d-9*dble(ao_num)**4 < dble(qp_max_mem)) then + call four_idx_dgemm + else + call add_integrals_to_map(full_ijkl_bitmask_4) + endif endif call wall_time(wall_2) @@ -77,6 +81,94 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] END_PROVIDER +subroutine four_idx_dgemm + implicit none + integer :: p,q,r,s,i,j,k,l + double precision, allocatable :: a1(:,:,:,:) + double precision, allocatable :: a2(:,:,:,:) + + allocate (a1(ao_num,ao_num,ao_num,ao_num)) + + print *, 'Getting AOs' + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s) + do s=1,ao_num + do r=1,ao_num + do q=1,ao_num + call get_ao_two_e_integrals(q,r,s,ao_num,a1(1,q,r,s)) + enddo + enddo + enddo + !$OMP END PARALLEL DO + + print *, '1st transformation' + ! 1st transformation + allocate (a2(ao_num,ao_num,ao_num,mo_num)) + call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num)) + + ! 2nd transformation + print *, '2nd transformation' + deallocate (a1) + allocate (a1(ao_num,ao_num,mo_num,mo_num)) + call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num)) + + ! 3rd transformation + print *, '3rd transformation' + deallocate (a2) + allocate (a2(ao_num,mo_num,mo_num,mo_num)) + call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num)) + + ! 4th transformation + print *, '4th transformation' + deallocate (a1) + allocate (a1(mo_num,mo_num,mo_num,mo_num)) + call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num)) + + deallocate (a2) + + integer :: n_integrals, size_buffer + integer(key_kind) , allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_value(:) + size_buffer = min(ao_num*ao_num*ao_num,16000000) + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals) + allocate ( buffer_i(size_buffer), buffer_value(size_buffer) ) + + n_integrals = 0 + !$OMP DO + do l=1,mo_num + do k=1,mo_num + do j=1,l + do i=1,k + if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = a1(i,j,k,l) + !DIR$ FORCEINLINE + call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + !$OMP END DO + + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + + deallocate (a1) + + call map_unique(mo_integrals_map) + + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + +end subroutine subroutine add_integrals_to_map(mask_ijkl) use bitmasks diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f index 6c66c8ec..64c15bf7 100644 --- a/src/tools/truncate_wf.irp.f +++ b/src/tools/truncate_wf.irp.f @@ -54,11 +54,23 @@ subroutine routine_s2 double precision, allocatable :: psi_coef_tmp(:,:) integer :: i,j,k double precision :: accu(N_states) + integer :: weights(0:16), ix + double precision :: x - print *, 'Weights of the CFG' + weights(:) = 0 do i=1,N_det - print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:))) + x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0) + ix = min(int(x), 16) + weights(ix) += 1 enddo + + print *, 'Histogram of the weights of the CFG' + do i=0,15 + print *, ' 10^{-', i, '} ', weights(i) + end do + print *, '< 10^{-', 15, '} ', weights(16) + + print*, 'Min weight of the configuration?' read(5,*) wmin diff --git a/src/utils/format_w_error.irp.f b/src/utils/format_w_error.irp.f new file mode 100644 index 00000000..1378d367 --- /dev/null +++ b/src/utils/format_w_error.irp.f @@ -0,0 +1,71 @@ +subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error) + + implicit none + + BEGIN_DOC + ! Format for double precision, value(error) + END_DOC + + ! in + ! | value | double precision | value... | + ! | error | double precision | error... | + ! | size_nb | integer | X in FX.Y | + ! | max_nb_digits | integer | Max Y in FX.Y | + + ! out + ! | format_value | character | string FX.Y for the format | + ! | str_error | character | string of the error | + + ! internal + ! | str_size | character | size in string format | + ! | nb_digits | integer | number of digits Y in FX.Y depending of the error | + ! | str_nb_digits | character | nb_digits in string format | + ! | str_exp | character | string of the value in exponential format | + + ! in + double precision, intent(in) :: error, value + integer, intent(in) :: size_nb, max_nb_digits + + ! out + character(len=20), intent(out) :: str_error, format_value + + ! internal + character(len=20) :: str_size, str_nb_digits, str_exp + integer :: nb_digits + + ! max_nb_digit: Y max + ! size_nb = Size of the double: X (FX.Y) + write(str_size,'(I3)') size_nb + + ! Error + write(str_exp,'(1pE20.0)') error + str_error = trim(adjustl(str_exp)) + + ! Number of digit: Y (FX.Y) from the exponent + str_nb_digits = str_exp(19:20) + read(str_nb_digits,*) nb_digits + + ! If the error is 0d0 + if (error <= 1d-16) then + write(str_nb_digits,*) max_nb_digits + endif + + ! If the error is too small + if (nb_digits > max_nb_digits) then + write(str_nb_digits,*) max_nb_digits + str_error(1:1) = '0' + endif + + ! If the error is too big (>= 0.5) + if (error >= 0.5d0) then + str_nb_digits = '1' + str_error(1:1) = '*' + endif + + ! FX.Y,A1,A1,A1 for value(str_error) + !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1' + + ! FX.Y just for the value + format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits)) + +end diff --git a/src/utils/map_module.f90 b/src/utils/map_module.f90 index 98e73470..ceaec874 100644 --- a/src/utils/map_module.f90 +++ b/src/utils/map_module.f90 @@ -238,11 +238,11 @@ subroutine cache_map_sort(map) iorder(i) = i enddo if (cache_key_kind == 2) then - call i2radix_sort(map%key,iorder,map%n_elements,-1) + call i2sort(map%key,iorder,map%n_elements,-1) else if (cache_key_kind == 4) then - call iradix_sort(map%key,iorder,map%n_elements,-1) + call isort(map%key,iorder,map%n_elements,-1) else if (cache_key_kind == 8) then - call i8radix_sort(map%key,iorder,map%n_elements,-1) + call i8sort(map%key,iorder,map%n_elements,-1) endif if (integral_kind == 4) then call set_order(map%value,iorder,map%n_elements) diff --git a/src/utils/qsort.c b/src/utils/qsort.c new file mode 100644 index 00000000..c011b35a --- /dev/null +++ b/src/utils/qsort.c @@ -0,0 +1,373 @@ +/* [[file:~/qp2/src/utils/qsort.org::*Generated%20C%20file][Generated C file:1]] */ +#include +#include + +struct int16_t_comp { + int16_t x; + int32_t i; +}; + +int compare_int16_t( const void * l, const void * r ) +{ + const int16_t * restrict _l= l; + const int16_t * restrict _r= r; + if( *_l > *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int16_t(int16_t* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct int16_t_comp* A = malloc(isize * sizeof(struct int16_t_comp)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int16_t_big(int16_t* restrict A_in, int64_t* restrict iorder, int64_t isize) { + struct int16_t_comp_big* A = malloc(isize * sizeof(struct int16_t_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int32_t(int32_t* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct int32_t_comp* A = malloc(isize * sizeof(struct int32_t_comp)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int32_t_big(int32_t* restrict A_in, int64_t* restrict iorder, int64_t isize) { + struct int32_t_comp_big* A = malloc(isize * sizeof(struct int32_t_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int64_t(int64_t* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct int64_t_comp* A = malloc(isize * sizeof(struct int64_t_comp)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_int64_t_big(int64_t* restrict A_in, int64_t* restrict iorder, int64_t isize) { + struct int64_t_comp_big* A = malloc(isize * sizeof(struct int64_t_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_double(double* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct double_comp* A = malloc(isize * sizeof(struct double_comp)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_double_big(double* restrict A_in, int64_t* restrict iorder, int64_t isize) { + struct double_comp_big* A = malloc(isize * sizeof(struct double_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_float(float* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct float_comp* A = malloc(isize * sizeof(struct float_comp)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_float_big(float* restrict A_in, int64_t* restrict iorder, int64_t isize) { + struct float_comp_big* A = malloc(isize * sizeof(struct float_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i *_r ) return 1; + if( *_l < *_r ) return -1; + return 0; +} + +void qsort_TYPE_big(TYPE* restrict A_in, int32_t* restrict iorder, int32_t isize) { + struct TYPE_comp_big* A = malloc(isize * sizeof(struct TYPE_comp_big)); + if (A == NULL) return; + + for (int i=0 ; i> +""" +for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]: + print( data.replace("TYPE", typ).replace("_big", "") ) + print( data.replace("int32_t", "int64_t").replace("TYPE", typ) ) +#+end_src + +#+NAME: replaced_f +#+begin_src python :results output :noweb yes +data = """ +<> +""" +c1 = { + "int16_t": "i2", + "int32_t": "i", + "int64_t": "i8", + "double": "d", + "float": "" +} +c2 = { + "int16_t": "integer", + "int32_t": "integer", + "int64_t": "integer", + "double": "real", + "float": "real" +} + +for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]: + print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") ) + print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) ) +#+end_src + +#+NAME: replaced_f2 +#+begin_src python :results output :noweb yes +data = """ +<> +""" +c1 = { + "int16_t": "i2", + "int32_t": "i", + "int64_t": "i8", + "double": "d", + "float": "" +} +c2 = { + "int16_t": "integer", + "int32_t": "integer", + "int64_t": "integer", + "double": "real", + "float": "real" +} + +for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]: + print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") ) + print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) ) +#+end_src + +* Generated C file + +#+BEGIN_SRC c :comments link :tangle qsort.c :noweb yes +#include +#include +<> +#+END_SRC + +* Generated Fortran file + +#+BEGIN_SRC f90 :tangle qsort_module.f90 :noweb yes +module qsort_module + use iso_c_binding + + interface + <> + end interface + +end module qsort_module + +<> + +#+END_SRC + diff --git a/src/utils/qsort_module.f90 b/src/utils/qsort_module.f90 new file mode 100644 index 00000000..a72a4f9e --- /dev/null +++ b/src/utils/qsort_module.f90 @@ -0,0 +1,347 @@ +module qsort_module + use iso_c_binding + + interface + + subroutine i2sort_c(A, iorder, isize) bind(C, name="qsort_int16_t") + use iso_c_binding + integer(c_int32_t), value :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int16_t) :: A(isize) + end subroutine i2sort_c + + subroutine i2sort_noidx_c(A, isize) bind(C, name="qsort_int16_t_noidx") + use iso_c_binding + integer(c_int32_t), value :: isize + integer (c_int16_t) :: A(isize) + end subroutine i2sort_noidx_c + + + + subroutine i2sort_big_c(A, iorder, isize) bind(C, name="qsort_int16_t_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int16_t) :: A(isize) + end subroutine i2sort_big_c + + subroutine i2sort_noidx_big_c(A, isize) bind(C, name="qsort_int16_t_noidx_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer (c_int16_t) :: A(isize) + end subroutine i2sort_noidx_big_c + + + + subroutine isort_c(A, iorder, isize) bind(C, name="qsort_int32_t") + use iso_c_binding + integer(c_int32_t), value :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int32_t) :: A(isize) + end subroutine isort_c + + subroutine isort_noidx_c(A, isize) bind(C, name="qsort_int32_t_noidx") + use iso_c_binding + integer(c_int32_t), value :: isize + integer (c_int32_t) :: A(isize) + end subroutine isort_noidx_c + + + + subroutine isort_big_c(A, iorder, isize) bind(C, name="qsort_int32_t_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int32_t) :: A(isize) + end subroutine isort_big_c + + subroutine isort_noidx_big_c(A, isize) bind(C, name="qsort_int32_t_noidx_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer (c_int32_t) :: A(isize) + end subroutine isort_noidx_big_c + + + + subroutine i8sort_c(A, iorder, isize) bind(C, name="qsort_int64_t") + use iso_c_binding + integer(c_int32_t), value :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int64_t) :: A(isize) + end subroutine i8sort_c + + subroutine i8sort_noidx_c(A, isize) bind(C, name="qsort_int64_t_noidx") + use iso_c_binding + integer(c_int32_t), value :: isize + integer (c_int64_t) :: A(isize) + end subroutine i8sort_noidx_c + + + + subroutine i8sort_big_c(A, iorder, isize) bind(C, name="qsort_int64_t_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int64_t) :: A(isize) + end subroutine i8sort_big_c + + subroutine i8sort_noidx_big_c(A, isize) bind(C, name="qsort_int64_t_noidx_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer (c_int64_t) :: A(isize) + end subroutine i8sort_noidx_big_c + + + + subroutine dsort_c(A, iorder, isize) bind(C, name="qsort_double") + use iso_c_binding + integer(c_int32_t), value :: isize + integer(c_int32_t) :: iorder(isize) + real (c_double) :: A(isize) + end subroutine dsort_c + + subroutine dsort_noidx_c(A, isize) bind(C, name="qsort_double_noidx") + use iso_c_binding + integer(c_int32_t), value :: isize + real (c_double) :: A(isize) + end subroutine dsort_noidx_c + + + + subroutine dsort_big_c(A, iorder, isize) bind(C, name="qsort_double_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer(c_int64_t) :: iorder(isize) + real (c_double) :: A(isize) + end subroutine dsort_big_c + + subroutine dsort_noidx_big_c(A, isize) bind(C, name="qsort_double_noidx_big") + use iso_c_binding + integer(c_int64_t), value :: isize + real (c_double) :: A(isize) + end subroutine dsort_noidx_big_c + + + + subroutine sort_c(A, iorder, isize) bind(C, name="qsort_float") + use iso_c_binding + integer(c_int32_t), value :: isize + integer(c_int32_t) :: iorder(isize) + real (c_float) :: A(isize) + end subroutine sort_c + + subroutine sort_noidx_c(A, isize) bind(C, name="qsort_float_noidx") + use iso_c_binding + integer(c_int32_t), value :: isize + real (c_float) :: A(isize) + end subroutine sort_noidx_c + + + + subroutine sort_big_c(A, iorder, isize) bind(C, name="qsort_float_big") + use iso_c_binding + integer(c_int64_t), value :: isize + integer(c_int64_t) :: iorder(isize) + real (c_float) :: A(isize) + end subroutine sort_big_c + + subroutine sort_noidx_big_c(A, isize) bind(C, name="qsort_float_noidx_big") + use iso_c_binding + integer(c_int64_t), value :: isize + real (c_float) :: A(isize) + end subroutine sort_noidx_big_c + + + + end interface + +end module qsort_module + + +subroutine i2sort(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int32_t) :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int16_t) :: A(isize) + call i2sort_c(A, iorder, isize) +end subroutine i2sort + +subroutine i2sort_noidx(A, isize) + use iso_c_binding + use qsort_module + integer(c_int32_t) :: isize + integer (c_int16_t) :: A(isize) + call i2sort_noidx_c(A, isize) +end subroutine i2sort_noidx + + + +subroutine i2sort_big(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int64_t) :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int16_t) :: A(isize) + call i2sort_big_c(A, iorder, isize) +end subroutine i2sort_big + +subroutine i2sort_noidx_big(A, isize) + use iso_c_binding + use qsort_module + integer(c_int64_t) :: isize + integer (c_int16_t) :: A(isize) + call i2sort_noidx_big_c(A, isize) +end subroutine i2sort_noidx_big + + + +subroutine isort(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int32_t) :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int32_t) :: A(isize) + call isort_c(A, iorder, isize) +end subroutine isort + +subroutine isort_noidx(A, isize) + use iso_c_binding + use qsort_module + integer(c_int32_t) :: isize + integer (c_int32_t) :: A(isize) + call isort_noidx_c(A, isize) +end subroutine isort_noidx + + + +subroutine isort_big(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int64_t) :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int32_t) :: A(isize) + call isort_big_c(A, iorder, isize) +end subroutine isort_big + +subroutine isort_noidx_big(A, isize) + use iso_c_binding + use qsort_module + integer(c_int64_t) :: isize + integer (c_int32_t) :: A(isize) + call isort_noidx_big_c(A, isize) +end subroutine isort_noidx_big + + + +subroutine i8sort(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int32_t) :: isize + integer(c_int32_t) :: iorder(isize) + integer (c_int64_t) :: A(isize) + call i8sort_c(A, iorder, isize) +end subroutine i8sort + +subroutine i8sort_noidx(A, isize) + use iso_c_binding + use qsort_module + integer(c_int32_t) :: isize + integer (c_int64_t) :: A(isize) + call i8sort_noidx_c(A, isize) +end subroutine i8sort_noidx + + + +subroutine i8sort_big(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int64_t) :: isize + integer(c_int64_t) :: iorder(isize) + integer (c_int64_t) :: A(isize) + call i8sort_big_c(A, iorder, isize) +end subroutine i8sort_big + +subroutine i8sort_noidx_big(A, isize) + use iso_c_binding + use qsort_module + integer(c_int64_t) :: isize + integer (c_int64_t) :: A(isize) + call i8sort_noidx_big_c(A, isize) +end subroutine i8sort_noidx_big + + + +subroutine dsort(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int32_t) :: isize + integer(c_int32_t) :: iorder(isize) + real (c_double) :: A(isize) + call dsort_c(A, iorder, isize) +end subroutine dsort + +subroutine dsort_noidx(A, isize) + use iso_c_binding + use qsort_module + integer(c_int32_t) :: isize + real (c_double) :: A(isize) + call dsort_noidx_c(A, isize) +end subroutine dsort_noidx + + + +subroutine dsort_big(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int64_t) :: isize + integer(c_int64_t) :: iorder(isize) + real (c_double) :: A(isize) + call dsort_big_c(A, iorder, isize) +end subroutine dsort_big + +subroutine dsort_noidx_big(A, isize) + use iso_c_binding + use qsort_module + integer(c_int64_t) :: isize + real (c_double) :: A(isize) + call dsort_noidx_big_c(A, isize) +end subroutine dsort_noidx_big + + + +subroutine sort(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int32_t) :: isize + integer(c_int32_t) :: iorder(isize) + real (c_float) :: A(isize) + call sort_c(A, iorder, isize) +end subroutine sort + +subroutine sort_noidx(A, isize) + use iso_c_binding + use qsort_module + integer(c_int32_t) :: isize + real (c_float) :: A(isize) + call sort_noidx_c(A, isize) +end subroutine sort_noidx + + + +subroutine sort_big(A, iorder, isize) + use qsort_module + use iso_c_binding + integer(c_int64_t) :: isize + integer(c_int64_t) :: iorder(isize) + real (c_float) :: A(isize) + call sort_big_c(A, iorder, isize) +end subroutine sort_big + +subroutine sort_noidx_big(A, isize) + use iso_c_binding + use qsort_module + integer(c_int64_t) :: isize + real (c_float) :: A(isize) + call sort_noidx_big_c(A, isize) +end subroutine sort_noidx_big diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index ff40263c..089c3871 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -1,222 +1,4 @@ BEGIN_TEMPLATE - subroutine insertion_$Xsort (x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize) using the insertion sort algorithm. - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - $type :: xtmp - integer :: i, i0, j, jmax - - do i=2,isize - xtmp = x(i) - i0 = iorder(i) - j=i-1 - do while (j>0) - if ((x(j) <= xtmp)) exit - x(j+1) = x(j) - iorder(j+1) = iorder(j) - j=j-1 - enddo - x(j+1) = xtmp - iorder(j+1) = i0 - enddo - end subroutine insertion_$Xsort - - subroutine quick_$Xsort(x, iorder, isize) - implicit none - BEGIN_DOC - ! Sort array x(isize) using the quicksort algorithm. - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer, external :: omp_get_num_threads - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - end - - recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) - implicit none - integer, intent(in) :: isize, first, last, level - integer,intent(inout) :: iorder(isize) - $type, intent(inout) :: x(isize) - $type :: c, tmp - integer :: itmp - integer :: i, j - - if(isize<2)return - - c = x( shiftr(first+last,1) ) - i = first - j = last - do - do while (x(i) < c) - i=i+1 - end do - do while (c < x(j)) - j=j-1 - end do - if (i >= j) exit - tmp = x(i) - x(i) = x(j) - x(j) = tmp - itmp = iorder(i) - iorder(i) = iorder(j) - iorder(j) = itmp - i=i+1 - j=j-1 - enddo - if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then - if (first < i-1) then - call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) - endif - if (j+1 < last) then - call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) - endif - else - if (first < i-1) then - call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) - endif - if (j+1 < last) then - call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) - endif - endif - end - - subroutine heap_$Xsort(x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize) using the heap sort algorithm. - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - - integer :: i, k, j, l, i0 - $type :: xtemp - - l = isize/2+1 - k = isize - do while (.True.) - if (l>1) then - l=l-1 - xtemp = x(l) - i0 = iorder(l) - else - xtemp = x(k) - i0 = iorder(k) - x(k) = x(1) - iorder(k) = iorder(1) - k = k-1 - if (k == 1) then - x(1) = xtemp - iorder(1) = i0 - exit - endif - endif - i=l - j = shiftl(l,1) - do while (j1) then - l=l-1 - xtemp = x(l) - i0 = iorder(l) - else - xtemp = x(k) - i0 = iorder(k) - x(k) = x(1) - iorder(k) = iorder(1) - k = k-1 - if (k == 1) then - x(1) = xtemp - iorder(1) = i0 - exit - endif - endif - i=l - j = shiftl(l,1) - do while (j0_8) - if (x(j)<=xtmp) exit - x(j+1_8) = x(j) - iorder(j+1_8) = iorder(j) - j = j-1_8 - enddo - x(j+1_8) = xtmp - iorder(j+1_8) = i0 - enddo - - end subroutine insertion_$Xsort_big - subroutine $Xset_order_big(x,iorder,isize) implicit none BEGIN_DOC @@ -565,223 +90,3 @@ SUBST [ X, type ] END_TEMPLATE -BEGIN_TEMPLATE - -recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) - implicit none - - BEGIN_DOC - ! Sort integer array x(isize) using the radix sort algorithm. - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - ! iradix should be -1 in input. - END_DOC - integer*$int_type, intent(in) :: isize - integer*$int_type, intent(inout) :: iorder(isize) - integer*$type, intent(inout) :: x(isize) - integer, intent(in) :: iradix - integer :: iradix_new - integer*$type, allocatable :: x2(:), x1(:) - integer*$type :: i4 ! data type - integer*$int_type, allocatable :: iorder1(:),iorder2(:) - integer*$int_type :: i0, i1, i2, i3, i ! index type - integer*$type :: mask - integer :: err - !DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1 - - if (isize < 2) then - return - endif - - if (iradix == -1) then ! Sort Positive and negative - - allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays' - stop - endif - - i1=1_$int_type - i2=1_$int_type - do i=1_$int_type,isize - if (x(i) < 0_$type) then - iorder1(i1) = iorder(i) - x1(i1) = -x(i) - i1 = i1+1_$int_type - else - iorder2(i2) = iorder(i) - x2(i2) = x(i) - i2 = i2+1_$int_type - endif - enddo - i1=i1-1_$int_type - i2=i2-1_$int_type - - do i=1_$int_type,i2 - iorder(i1+i) = iorder2(i) - x(i1+i) = x2(i) - enddo - deallocate(x2,iorder2,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x2, iorder2' - stop - endif - - - if (i1 > 1_$int_type) then - call $Xradix_sort$big(x1,iorder1,i1,-2) - do i=1_$int_type,i1 - x(i) = -x1(1_$int_type+i1-i) - iorder(i) = iorder1(1_$int_type+i1-i) - enddo - endif - - if (i2>1_$int_type) then - call $Xradix_sort$big(x(i1+1_$int_type),iorder(i1+1_$int_type),i2,-2) - endif - - deallocate(x1,iorder1,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x1, iorder1' - stop - endif - return - - else if (iradix == -2) then ! Positive - - ! Find most significant bit - - i0 = 0_$int_type - i4 = maxval(x) - - iradix_new = max($integer_size-1-leadz(i4),1) - mask = ibset(0_$type,iradix_new) - - allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays' - stop - endif - - i1=1_$int_type - i2=1_$int_type - - do i=1_$int_type,isize - if (iand(mask,x(i)) == 0_$type) then - iorder1(i1) = iorder(i) - x1(i1) = x(i) - i1 = i1+1_$int_type - else - iorder2(i2) = iorder(i) - x2(i2) = x(i) - i2 = i2+1_$int_type - endif - enddo - i1=i1-1_$int_type - i2=i2-1_$int_type - - do i=1_$int_type,i1 - iorder(i0+i) = iorder1(i) - x(i0+i) = x1(i) - enddo - i0 = i0+i1 - i3 = i0 - deallocate(x1,iorder1,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x1, iorder1' - stop - endif - - - do i=1_$int_type,i2 - iorder(i0+i) = iorder2(i) - x(i0+i) = x2(i) - enddo - i0 = i0+i2 - deallocate(x2,iorder2,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x2, iorder2' - stop - endif - - - if (i3>1_$int_type) then - call $Xradix_sort$big(x,iorder,i3,iradix_new-1) - endif - - if (isize-i3>1_$int_type) then - call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) - endif - - return - endif - - ASSERT (iradix >= 0) - - if (isize < 48) then - call insertion_$Xsort$big(x,iorder,isize) - return - endif - - - allocate(x2(isize),iorder2(isize),stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays x1, iorder1' - stop - endif - - - mask = ibset(0_$type,iradix) - i0=1_$int_type - i1=1_$int_type - - do i=1_$int_type,isize - if (iand(mask,x(i)) == 0_$type) then - iorder(i0) = iorder(i) - x(i0) = x(i) - i0 = i0+1_$int_type - else - iorder2(i1) = iorder(i) - x2(i1) = x(i) - i1 = i1+1_$int_type - endif - enddo - i0=i0-1_$int_type - i1=i1-1_$int_type - - do i=1_$int_type,i1 - iorder(i0+i) = iorder2(i) - x(i0+i) = x2(i) - enddo - - deallocate(x2,iorder2,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays x2, iorder2' - stop - endif - - - if (iradix == 0) then - return - endif - - - if (i1>1_$int_type) then - call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) - endif - if (i0>1) then - call $Xradix_sort$big(x,iorder,i0,iradix-1) - endif - - end - -SUBST [ X, type, integer_size, is_big, big, int_type ] - i ; 4 ; 32 ; .False. ; ; 4 ;; - i8 ; 8 ; 64 ; .False. ; ; 4 ;; - i2 ; 2 ; 16 ; .False. ; ; 4 ;; - i ; 4 ; 32 ; .True. ; _big ; 8 ;; - i8 ; 8 ; 64 ; .True. ; _big ; 8 ;; -END_TEMPLATE - - - diff --git a/src/utils/units.irp.f b/src/utils/units.irp.f new file mode 100644 index 00000000..1850b28b --- /dev/null +++ b/src/utils/units.irp.f @@ -0,0 +1,22 @@ +BEGIN_PROVIDER [double precision, ha_to_ev] + + implicit none + BEGIN_DOC + ! Converstion from Hartree to eV + END_DOC + + ha_to_ev = 27.211396641308d0 + +END_PROVIDER + +BEGIN_PROVIDER [double precision, au_to_D] + + implicit none + BEGIN_DOC + ! Converstion from au to Debye + END_DOC + + au_to_D = 2.5415802529d0 + +END_PROVIDER + diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index ef846bdb..84593031 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -37,6 +37,10 @@ double precision function binom_func(i,j) else binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) ) endif + + ! To avoid .999999 numbers + binom_func = floor(binom_func + 0.5d0) + end diff --git a/tests/bats/common.bats.sh b/tests/bats/common.bats.sh index f6ea4023..802c0232 100644 --- a/tests/bats/common.bats.sh +++ b/tests/bats/common.bats.sh @@ -46,7 +46,7 @@ function test_exe() { run_only_test() { if [[ "$BATS_TEST_DESCRIPTION" != "$1" ]] && [[ "$BATS_TEST_NUMBER" != "$1" ]]; then - if [[ -z $BATS_TEST_FILENAME ]] ; then + if [[ -z "$BATS_TEST_FILENAME" ]] ; then exit 0 else skip diff --git a/travis/compilation.sh b/travis/compilation.sh deleted file mode 100755 index 071b4872..00000000 --- a/travis/compilation.sh +++ /dev/null @@ -1,16 +0,0 @@ -#!/bin/bash -# Stage 2 - -# Extract cache from config stage -cd ../ -tar -zxf $HOME/cache/config.tgz - -# Configure QP2 -cd qp2 -source ./quantum_package.rc -ninja -j 1 -v || exit -1 - -# Create cache -cd .. -tar -zcf $HOME/cache/compil.tgz qp2 && rm $HOME/cache/config.tgz - diff --git a/travis/configuration.sh b/travis/configuration.sh deleted file mode 100755 index f925107d..00000000 --- a/travis/configuration.sh +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -# Stage 1 - -# Configure QP2 -./configure --download all --install all --config ./config/travis.cfg || exit -1 - -# Create cache -cd ../ -tar -zcf $HOME/cache/config.tgz qp2 - diff --git a/travis/testing.sh b/travis/testing.sh deleted file mode 100755 index f67bd106..00000000 --- a/travis/testing.sh +++ /dev/null @@ -1,16 +0,0 @@ -#!/bin/bash -# Stage 3 - -# Extract cache from compile stage -cd ../ -tar -zxf $HOME/cache/compil.tgz - -# Configure QP2 -cd qp2 -source ./quantum_package.rc -exec qp_test -a && rm $HOME/cache/compil.tgz - - - - -