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/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 091423e4..07ad2236 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -224,7 +224,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_set_frozen_core b/bin/qp_set_frozen_core index 879c71de..bc6f6834 100755 --- a/bin/qp_set_frozen_core +++ b/bin/qp_set_frozen_core @@ -7,12 +7,13 @@ setting all MOs as Active, except the n/2 first ones which are set as Core. If pseudo-potentials are used, all the MOs are set as Active. Usage: - qp_set_frozen_core [-q|--query] [(-l|-s|--large|--small)] EZFIO_DIR + qp_set_frozen_core [-q|--query] [(-l|-s|-u|--large|--small|--unset)] EZFIO_DIR Options: -q --query Prints in the standard output the number of frozen MOs -l --large Use a small core -s --small Use a large core + -u --unset Unset frozen core Default numbers of frozen electrons: @@ -88,7 +89,9 @@ def main(arguments): elif charge <= 54: n_frozen += 9 elif charge <= 86: n_frozen += 18 elif charge <= 118: n_frozen += 27 + elif arguments["--unset"]: + n_frozen = 0 else: # default for charge in ezfio.nuclei_nucl_charge: if charge <= 4: pass 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/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_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/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 33a3bc89..336161f8 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -1,7 +1,7 @@ ! Spherical to cartesian transformation matrix obtained with ! Horton (http://theochem.github.com/horton/, 2015) -! First index is the index of the carteisan AO, obtained by ao_power_index +! First index is the index of the cartesian AO, obtained by ao_power_index ! Second index is the index of the spherical AO BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] 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 14b1d060..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) @@ -523,28 +541,74 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) if(c > 2) then eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % pt2(pt2_stoch_istate) = eqt eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % variance(pt2_stoch_istate) = eqt eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0)) + eqta(:) = dsqrt(eqta(:) / (dble(c) - 1.5d0)) pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) 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.4)', 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 @@ -842,9 +906,8 @@ END_PROVIDER do t=1, pt2_N_teeth tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))) endif - ASSERT(tooth_width > 0.d0) do i=pt2_n_0(t)+1, pt2_n_0(t+1) pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width end do diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index f1001f89..30fc7ce0 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -116,10 +116,10 @@ subroutine run_pt2_slave_small(thread,iproc,energy) do k=1,n_tasks call pt2_alloc(pt2_data(k),N_states) b%cur = 0 - double precision :: time2 - call wall_time(time2) +! double precision :: time2 +! call wall_time(time2) call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) - call wall_time(time1) +! call wall_time(time1) ! print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) enddo call wall_time(time1) @@ -190,8 +190,12 @@ subroutine run_pt2_slave_large(thread,iproc,energy) integer :: bsize ! Size of selection buffers logical :: sending + double precision :: time_shift + PROVIDE global_selection_buffer global_selection_buffer_lock + call random_number(time_shift) + time_shift = time_shift*15.d0 zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() @@ -209,6 +213,9 @@ subroutine run_pt2_slave_large(thread,iproc,energy) sending = .False. done = .False. + double precision :: time0, time1 + call wall_time(time0) + time0 = time0+time_shift do while (.not.done) integer, external :: get_tasks_from_taskserver @@ -244,19 +251,24 @@ subroutine run_pt2_slave_large(thread,iproc,energy) done = .true. endif call sort_selection_buffer(b) - call omp_set_lock(global_selection_buffer_lock) - global_selection_buffer%mini = b%mini - call merge_selection_buffers(b,global_selection_buffer) - b%cur=0 - call omp_unset_lock(global_selection_buffer_lock) + + call wall_time(time1) +! if (time1-time0 > 15.d0) then + call omp_set_lock(global_selection_buffer_lock) + global_selection_buffer%mini = b%mini + call merge_selection_buffers(b,global_selection_buffer) + b%cur=0 + call omp_unset_lock(global_selection_buffer_lock) + call wall_time(time0) +! endif + + call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) if ( iproc == 1 .or. i_generator < 100 .or. done) then call omp_set_lock(global_selection_buffer_lock) - call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending) global_selection_buffer%cur = 0 call omp_unset_lock(global_selection_buffer_lock) else - call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), b, (/task_id/), 1,sending) endif 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 f1ec6ff6..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)) @@ -474,17 +477,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d ! endif do i=1,fullinteresting(0) - do k=1,N_int - fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) - fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) - enddo + fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i)) enddo do i=1,interesting(0) - do k=1,N_int - minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) - minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) - enddo + minilist(:,:,i) = psi_det_sorted(:,:,interesting(i)) enddo do s2=s1,2 @@ -572,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 @@ -713,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 @@ -783,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) @@ -834,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 @@ -1556,7 +1594,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi/selection_buffer.irp.f index 79899139..1f743e0e 100644 --- a/src/cipsi/selection_buffer.irp.f +++ b/src/cipsi/selection_buffer.irp.f @@ -92,38 +92,51 @@ subroutine merge_selection_buffers(b1, b2) allocate(val(sze), detmp(N_int, 2, sze)) i1=1 i2=1 - do i=1,nmwen - if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then - exit - else if (i1 > b1%cur) then - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 - else if (i2 > b2%cur) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 - else - if (b1%val(i1) <= b2%val(i2)) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 + + select case (N_int) +BEGIN_TEMPLATE + case $case + do i=1,nmwen + if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then + exit + else if (i1 > b1%cur) then + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + else if (i2 > b2%cur) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 else - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 + if (b1%val(i1) <= b2%val(i2)) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 + else + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + endif endif - endif - enddo + enddo + do i=nmwen+1,b2%N + val(i) = 0.d0 +! detmp(1:$N_int,1,i) = 0_bit_kind +! detmp(1:$N_int,2,i) = 0_bit_kind + enddo +SUBST [ case, N_int ] +(1); 1;; +(2); 2;; +(3); 3;; +(4); 4;; +default; N_int;; +END_TEMPLATE + end select deallocate(b2%det, b2%val) - do i=nmwen+1,b2%N - val(i) = 0.d0 - detmp(1:N_int,1:2,i) = 0_bit_kind - enddo b2%det => detmp b2%val => val b2%mini = min(b2%mini,b2%val(b2%N)) diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index ab2294ad..2b16a5f7 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -62,6 +62,7 @@ subroutine run else call H_apply_cis endif + print*,'' print *, 'N_det = ', N_det print*,'******************************' print *, 'Energies of the states:' @@ -69,11 +70,13 @@ subroutine run print *, i, CI_energy(i) enddo if (N_states > 1) then - print*,'******************************' - print*,'Excitation energies ' + print*,'' + print*,'******************************************************' + print*,'Excitation energies (au) (eV)' do i = 2, N_states - print*, i ,CI_energy(i) - CI_energy(1) + print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev enddo + print*,'' endif call ezfio_set_cis_energy(CI_energy) 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..5f167686 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -69,7 +69,9 @@ subroutine run 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 +80,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 0e4dd91f..5190785a 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -341,6 +341,7 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl Get BFtoDeterminant Matrix ************************************/ + convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); int rowsI = 0; @@ -1331,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); @@ -1428,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); @@ -1696,7 +1700,6 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in gramSchmidt(overlapMatrixJ, rowsJ, colsJ, orthoMatrixJ); - int rowsA = 0; int colsA = 0; diff --git a/src/csf/configurations.irp.f b/src/csf/configurations.irp.f index c11a49a4..aebf53d9 100644 --- a/src/csf/configurations.irp.f +++ b/src/csf/configurations.irp.f @@ -856,6 +856,7 @@ end subroutine !end subroutine ! BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ] &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] implicit none @@ -944,6 +945,29 @@ end subroutine enddo deallocate(dets, old_order) + integer :: ndet_conf + do i = 1, N_configuration + ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1 + psi_configuration_n_det(i) = ndet_conf + enddo END_PROVIDER + +BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)] + implicit none + integer :: i,j,k,l + integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int) + n_elec_alpha_for_psi_configuration = 0 + do i = 1, N_configuration + j = psi_configuration_to_psi_det(2,i) + det_tmp(:,:) = psi_det(:,:,j) + k = 0 + do l = 1, N_int + det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i)) + k += popcnt(det_alpha(l)) + enddo + n_elec_alpha_for_psi_configuration(i) = k + enddo + +END_PROVIDER diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 9f7f4442..75f6e539 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -58,15 +58,14 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) enddo enddo - s = 0 + s = 0 ! s == total number of SOMOs do k=1,N_int if (psi_configuration(k,1,i) == 0_bit_kind) cycle s = s + popcnt(psi_configuration(k,1,i)) enddo - ! Test 1 - if(iand(MS,1) .EQ. 0) then - bfIcfg = max(1,nint((binom(s,s/2)-binom(s,s/2+1)))) + if(iand(s,1) .EQ. 0) then + bfIcfg = max(1,nint((binom(s,s/2)-binom(s,(s/2)+1)))) else bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) endif diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 6748b312..60865320 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -136,6 +136,7 @@ ncfgprev = cfg_seniority_index(i+2) end do print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration + END_PROVIDER diff --git a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f index 2f3d7f80..243e9995 100644 --- a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f +++ b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f @@ -1,5 +1,5 @@ -subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) +subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) use mmap_module implicit none BEGIN_DOC @@ -412,36 +412,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_d FREE nthreads_davidson end -subroutine hcalc_template(v,u,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! Template of routine for the application of H - ! - ! Here, it is done with the Hamiltonian matrix - ! - ! on the set of determinants of psi_det - ! - ! Computes $v = H | u \rangle$ - ! - END_DOC - integer, intent(in) :: N_st,sze - double precision, intent(in) :: u(sze,N_st) - double precision, intent(inout) :: v(sze,N_st) - integer :: i,j,istate - v = 0.d0 - do istate = 1, N_st - do i = 1, sze - do j = 1, sze - v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate) - enddo - enddo - do i = 1, sze - v(i,istate) += u(i,istate) * nuclear_repulsion - enddo - enddo -end - subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze) implicit none BEGIN_DOC 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/davidson/u0_hs2_u0.irp.f b/src/davidson/u0_hs2_u0.irp.f index 8f7bf06b..38fb56bd 100644 --- a/src/davidson/u0_hs2_u0.irp.f +++ b/src/davidson/u0_hs2_u0.irp.f @@ -203,7 +203,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, integer, allocatable :: doubles(:) integer, allocatable :: singles_a(:) integer, allocatable :: singles_b(:) - integer, allocatable :: idx(:), idx0(:) + integer, allocatable :: idx(:), buffer_lrow(:), idx0(:) integer :: maxab, n_singles_a, n_singles_b, kcol_prev integer*8 :: k8 logical :: compute_singles @@ -253,7 +253,7 @@ compute_singles=.True. !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & !$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, & !$OMP buffer, doubles, n_doubles, umax, & - !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & + !$OMP tmp_det2, hij, sij, idx, buffer_lrow, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, ratio, & !$OMP n_singles_b, k8, last_found,left,right,right_max) @@ -264,7 +264,7 @@ compute_singles=.True. singles_a(maxab), & singles_b(maxab), & doubles(maxab), & - idx(maxab), utl(N_st,block_size)) + idx(maxab), buffer_lrow(maxab), utl(N_st,block_size)) kcol_prev=-1 @@ -332,18 +332,20 @@ compute_singles=.True. l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) - !DIR$ UNROLL(8) - !DIR$ LOOP COUNT avg(50000) do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) ! hot spot + buffer_lrow(j) = lrow ASSERT (l_a <= N_det) idx(j) = l_a l_a = l_a+1 enddo + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, buffer_lrow(j)) ! hot spot + enddo j = j-1 call get_all_spin_singles_$N_int( & @@ -789,7 +791,7 @@ compute_singles=.True. end do !$OMP END DO - deallocate(buffer, singles_a, singles_b, doubles, idx, utl) + deallocate(buffer, singles_a, singles_b, doubles, idx, buffer_lrow, utl) !$OMP END PARALLEL end diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index b8c8658f..4b317025 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 = 1 + 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 = 1 + 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 8a5f1a2d..b411dda4 100644 --- a/src/determinants/dipole_moments.irp.f +++ b/src/determinants/dipole_moments.irp.f @@ -9,7 +9,7 @@ double precision :: weight, r(3) double precision :: cpu0,cpu1,nuclei_part_z,nuclei_part_y,nuclei_part_x - call cpu_time(cpu0) +! call cpu_time(cpu0) z_dipole_moment = 0.d0 y_dipole_moment = 0.d0 x_dipole_moment = 0.d0 @@ -26,10 +26,10 @@ enddo enddo - print*,'electron part for z_dipole = ',z_dipole_moment - print*,'electron part for y_dipole = ',y_dipole_moment - print*,'electron part for x_dipole = ',x_dipole_moment - +! print*,'electron part for z_dipole = ',z_dipole_moment +! print*,'electron part for y_dipole = ',y_dipole_moment +! print*,'electron part for x_dipole = ',x_dipole_moment +! nuclei_part_z = 0.d0 nuclei_part_y = 0.d0 nuclei_part_x = 0.d0 @@ -38,28 +38,43 @@ nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) enddo - print*,'nuclei part for z_dipole = ',nuclei_part_z - print*,'nuclei part for y_dipole = ',nuclei_part_y - print*,'nuclei part for x_dipole = ',nuclei_part_x - +! print*,'nuclei part for z_dipole = ',nuclei_part_z +! print*,'nuclei part for y_dipole = ',nuclei_part_y +! print*,'nuclei part for x_dipole = ',nuclei_part_x +! do istate = 1, N_states z_dipole_moment(istate) += nuclei_part_z y_dipole_moment(istate) += nuclei_part_y x_dipole_moment(istate) += nuclei_part_x enddo - call cpu_time(cpu1) - print*,'Time to provide the dipole moment :',cpu1-cpu0 +! call cpu_time(cpu1) +! print*,'Time to provide the dipole moment :',cpu1-cpu0 END_PROVIDER - subroutine print_z_dipole_moment_only + subroutine print_dipole_moments implicit none + integer :: i print*, '' print*, '' print*, '****************************************' - print*, 'z_dipole_moment = ',z_dipole_moment + write(*,'(A10)',advance='no') ' State : ' + do i = 1,N_states + write(*,'(i16)',advance='no') i + end do + write(*,*) '' + 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 print*, '****************************************' end diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f index d73b2dbf..2c1a8757 100644 --- a/src/determinants/s2.irp.f +++ b/src/determinants/s2.irp.f @@ -103,13 +103,17 @@ BEGIN_PROVIDER [ double precision, expected_s2] END_PROVIDER -BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] + BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] +&BEGIN_PROVIDER [ double precision, s_values, (N_states) ] implicit none BEGIN_DOC ! array of the averaged values of the S^2 operator on the various states END_DOC integer :: i call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size) + do i = 1, N_states + s_values(i) = 0.5d0 *(-1.d0 + dsqrt(1.d0 + 4 * s2_values(i))) + enddo END_PROVIDER 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/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 6fa6a4c7..38478d21 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -91,7 +91,19 @@ enddo END_PROVIDER - BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)] + BEGIN_PROVIDER [double precision, aos_lapl_in_r_array_transp, (ao_num, n_points_final_grid,3)] + implicit none + integer :: i,j,m + do i = 1, n_points_final_grid + do j = 1, ao_num + do m = 1, 3 + aos_lapl_in_r_array_transp(j,i,m) = aos_lapl_in_r_array(m,j,i) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)] implicit none BEGIN_DOC ! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point @@ -100,20 +112,20 @@ END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(ao_num,3) - double precision :: aos_lapl_array(ao_num,3) + double precision :: aos_grad_array(3,ao_num) + double precision :: aos_lapl_array(3,ao_num) !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points) - do m = 1, 3 - do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) - do j = 1, ao_num - aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) + do j = 1, ao_num + do m = 1, 3 + aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j) enddo enddo enddo diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 0a8b4d52..192cb25a 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -138,7 +138,7 @@ integer :: m mos_lapl_in_r_array = 0.d0 do m=1,3 - call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num) + call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array_transp(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num) enddo END_PROVIDER diff --git a/src/dressing/alpha_factory.irp.f b/src/dressing/alpha_factory.irp.f index 5eeeb1a6..c7adffe3 100644 --- a/src/dressing/alpha_factory.irp.f +++ b/src/dressing/alpha_factory.irp.f @@ -1179,7 +1179,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/ezfio_files/output.irp.f b/src/ezfio_files/output.irp.f index 48512f92..7b2663a0 100644 --- a/src/ezfio_files/output.irp.f +++ b/src/ezfio_files/output.irp.f @@ -25,7 +25,7 @@ subroutine write_time(iunit) ct = ct - output_cpu_time_0 call wall_time(wt) wt = wt - output_wall_time_0 - write(6,'(A,F14.6,A,F14.6,A)') & + write(6,'(A,F14.2,A,F14.2,A)') & '.. >>>>> [ WALL TIME: ', wt, ' s ] [ CPU TIME: ', ct, ' s ] <<<<< ..' write(6,*) end 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_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index cbf23a9a..fcbdde49 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -7,7 +7,7 @@ subroutine hcore_guess label = 'Guess' call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, & size(mo_one_e_integrals,1), & - size(mo_one_e_integrals,2),label,1,.false.) + size(mo_one_e_integrals,2),label,1,.true.) call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 ) call save_mos TOUCH mo_coef mo_label 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/print_dipole.irp.f b/src/tools/print_dipole.irp.f index 8351308e..8db9aa09 100644 --- a/src/tools/print_dipole.irp.f +++ b/src/tools/print_dipole.irp.f @@ -1,5 +1,7 @@ program print_dipole implicit none - call print_z_dipole_moment_only + read_wf = .True. + TOUCH read_wf + call print_dipole_moments end diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f index 7e51caaf..64eb1a1f 100644 --- a/src/tools/print_wf.irp.f +++ b/src/tools/print_wf.irp.f @@ -32,8 +32,9 @@ subroutine routine double precision :: norm_mono_a,norm_mono_b double precision :: norm_mono_a_2,norm_mono_b_2 double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2 - double precision :: norm_mono_a_pert,norm_mono_b_pert + double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1 double precision :: delta_e,coef_2_2 + norm_mono_a = 0.d0 norm_mono_b = 0.d0 norm_mono_a_2 = 0.d0 @@ -42,6 +43,7 @@ subroutine routine norm_mono_b_pert = 0.d0 norm_mono_a_pert_2 = 0.d0 norm_mono_b_pert_2 = 0.d0 + norm_double_1 = 0.d0 do i = 1, min(N_det_print_wf,N_det) print*,'' print*,'i = ',i @@ -93,6 +95,7 @@ subroutine routine print*,'h1,p1 = ',h1,p1 print*,'s2',s2 print*,'h2,p2 = ',h2,p2 + norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1)) endif print*,' = ',hij @@ -109,6 +112,7 @@ subroutine routine print*,'' print*,'L1 norm of mono alpha = ',norm_mono_a print*,'L1 norm of mono beta = ',norm_mono_b + print*,'L1 norm of double exc = ',norm_double_1 print*, '---' print*,'L2 norm of mono alpha = ',norm_mono_a_2 print*,'L2 norm of mono beta = ',norm_mono_b_2 diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index 4dadd2e6..a4dea15f 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -1,9 +1,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] implicit none BEGIN_DOC - ! two_e_dm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons - ! - ! + ! \sum_{\sigma \sigma'} + ! ! ! where the indices (i,j,k,l) belong to all MOs. ! @@ -12,7 +11,7 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO are set to zero ! The state-averaged two-electron energy : ! - ! \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < ii jj | kk ll > + ! \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < kk ll | ii jj > END_DOC two_e_dm_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb,istate 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/memory.irp.f b/src/utils/memory.irp.f index 3ea242b0..d5a066a1 100644 --- a/src/utils/memory.irp.f +++ b/src/utils/memory.irp.f @@ -114,7 +114,7 @@ subroutine print_memory_usage() call resident_memory(rss) call total_memory(mem) - write(*,'(A,F14.6,A,F14.6,A)') & + write(*,'(A,F14.3,A,F14.3,A)') & '.. >>>>> [ RES MEM : ', rss , & ' GB ] [ VIRT MEM : ', mem, ' GB ] <<<<< ..' end 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/set_multiple_levels_omp.irp.f b/src/utils/set_multiple_levels_omp.irp.f index b4764e4a..572a13f4 100644 --- a/src/utils/set_multiple_levels_omp.irp.f +++ b/src/utils/set_multiple_levels_omp.irp.f @@ -8,7 +8,7 @@ subroutine set_multiple_levels_omp(activate) logical, intent(in) :: activate if (activate) then - call omp_set_max_active_levels(5) + call omp_set_max_active_levels(3) IRP_IF SET_NESTED call omp_set_nested(.True.) 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 - - - - -