PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-15 10:33:50 +01:00
Code Releases Activity

Merge branch 'dev' into csf_verified

Browse Source
This commit is contained in:
vijay 2022-06-13 16:27:30 +02:00 committed by GitHub
commit 872ed7f283
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
77 changed files with 1891 additions and 1110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
.drone.yml Normal file
View File

@ -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]

2
bin/qp_convert_output_to_ezfio
View File

@ -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 #

5
bin/qp_set_frozen_core
View File

@ -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

9
bin/qp_test
View File

@ -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)

2
config/ifort_2019_avx.cfg
View File

@ -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

2
config/ifort_2019_avx_mpi.cfg
View File

@ -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

2
config/ifort_2019_mpi_rome.cfg
View File

@ -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

2
config/ifort_2019_rome.cfg
View File

@ -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

2
config/ifort_2019_sse4.cfg
View File

@ -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

2
config/ifort_2019_sse4_mpi.cfg
View File

@ -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

2
config/ifort_2019_xHost.cfg
View File

@ -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

2
config/ifort_2021_avx.cfg
View File

@ -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

2
config/ifort_2021_avx_mpi.cfg
View File

@ -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

2
config/ifort_2021_mpi_rome.cfg
View File

@ -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

2
config/ifort_2021_rome.cfg
View File

@ -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

2
config/ifort_2021_sse4.cfg
View File

@ -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

2
config/ifort_2021_sse4_mpi.cfg
View File

@ -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

2
config/ifort_2021_xHost.cfg
View File

@ -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

4
configure vendored
View File

@ -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

4
ocaml/Qputils.ml
View File

@ -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 )

139
ocaml/qp_tunnel.ml
View File

@ -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 <host> &
%!"
%!"
(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"

27
scripts/cipsi_save.sh Normal file
View File

@ -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

2
src/ao_basis/spherical_to_cartesian.irp.f
View File

@ -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) ]

6
src/cipsi/EZFIO.cfg
View File

@ -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

97
src/cipsi/pt2_stoch_routines.irp.f
View File

@ -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

32
src/cipsi/run_pt2_slave.irp.f
View File

@ -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

8
src/cipsi/run_selection_slave.irp.f
View File

@ -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

72
src/cipsi/selection.irp.f
View File

@ -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)

71
src/cipsi/selection_buffer.irp.f
View File

@ -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))

9
src/cis/cis.irp.f
View File

@ -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)

7
src/cis_read/EZFIO.cfg Normal file
View File

@ -0,0 +1,7 @@
[energy]
type: double precision
doc: Variational |CIS| energy
interface: ezfio
size: (determinants.n_states)

3
src/cis_read/NEED Normal file
View File

@ -0,0 +1,3 @@
selectors_full
generators_full
davidson_undressed

5
src/cis_read/README.rst Normal file
View File

@ -0,0 +1,5 @@
===
cis_read
===
Reads the input WF and performs all singles on top of it.

88
src/cis_read/cis_read.irp.f Normal file
View File

@ -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

14
src/cis_read/h_apply.irp.f Normal file
View File

@ -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

2
src/cisd/30.cisd.bats
View File

@ -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

59
src/cisd/cisd.irp.f
View File

@ -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

11
src/csf/cfgCI_utils.c
View File

@ -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;

24
src/csf/configurations.irp.f
View File

@ -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

7
src/csf/conversion.irp.f
View File

@ -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

1
src/csf/sigma_vector.irp.f
View File

@ -136,6 +136,7 @@
ncfgprev = cfg_seniority_index(i+2)
end do
print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
END_PROVIDER

32
src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
View File

@ -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

25
src/davidson/diagonalize_ci.irp.f
View File

@ -1,3 +1,13 @@
BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ]
implicit none
BEGIN_DOC
! If 'det', use <Psi_det|H|Psi_det> in Davidson
!
! If 'cfg', use <Psi_csf|H|Psi_csf> 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, &

16
src/davidson/u0_hs2_u0.irp.f
View File

@ -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

45
src/determinants/determinants.irp.f
View File

@ -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)

41
src/determinants/dipole_moments.irp.f
View File

@ -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

6
src/determinants/s2.irp.f
View File

@ -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

9
src/determinants/slater_rules.irp.f
View File

@ -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

3
src/determinants/slater_rules_wee_mono.irp.f
View File

@ -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)

34
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -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

2
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -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

2
src/dressing/alpha_factory.irp.f
View File

@ -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)

2
src/ezfio_files/output.irp.f
View File

@ -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

115
src/fci/40.fci.bats
View File

@ -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
}

6
src/fci/EZFIO.cfg
View File

@ -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)

3
src/iterations/print_extrapolation.irp.f
View File

@ -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

12
src/iterations/print_summary.irp.f
View File

@ -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

2
src/mo_guess/h_core_guess_routine.irp.f
View File

@ -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

12
src/mo_two_e_erf_ints/map_integrals_erf.irp.f
View File

@ -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)

94
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -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

4
src/tools/print_dipole.irp.f
View File

@ -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

6
src/tools/print_wf.irp.f
View File

@ -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*,'<Ref| H |D_I> = ',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

7
src/two_body_rdm/two_e_dm_mo.irp.f
View File

@ -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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
! \sum_{\sigma \sigma'}
! <Psi| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi>
!
! 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

71
src/utils/format_w_error.irp.f Normal file
View File

@ -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

6
src/utils/map_module.f90
View File

@ -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)

2
src/utils/memory.irp.f
View File

@ -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

373
src/utils/qsort.c Normal file
View File

@ -0,0 +1,373 @@
/* [[file:~/qp2/src/utils/qsort.org::*Generated%20C%20file][Generated C file:1]] */
#include <stdlib.h>
#include <stdint.h>
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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp), compare_int16_t);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int16_t_noidx(int16_t* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t);
}
struct int16_t_comp_big {
int16_t x;
int64_t i;
};
int compare_int16_t_big( 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_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp_big), compare_int16_t_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int16_t_noidx_big(int16_t* A, int64_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t_big);
}
struct int32_t_comp {
int32_t x;
int32_t i;
};
int compare_int32_t( const void * l, const void * r )
{
const int32_t * restrict _l= l;
const int32_t * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp), compare_int32_t);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int32_t_noidx(int32_t* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t);
}
struct int32_t_comp_big {
int32_t x;
int64_t i;
};
int compare_int32_t_big( const void * l, const void * r )
{
const int32_t * restrict _l= l;
const int32_t * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp_big), compare_int32_t_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int32_t_noidx_big(int32_t* A, int64_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t_big);
}
struct int64_t_comp {
int64_t x;
int32_t i;
};
int compare_int64_t( const void * l, const void * r )
{
const int64_t * restrict _l= l;
const int64_t * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp), compare_int64_t);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int64_t_noidx(int64_t* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t);
}
struct int64_t_comp_big {
int64_t x;
int64_t i;
};
int compare_int64_t_big( const void * l, const void * r )
{
const int64_t * restrict _l= l;
const int64_t * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp_big), compare_int64_t_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_int64_t_noidx_big(int64_t* A, int64_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t_big);
}
struct double_comp {
double x;
int32_t i;
};
int compare_double( const void * l, const void * r )
{
const double * restrict _l= l;
const double * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct double_comp), compare_double);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_double_noidx(double* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(double), compare_double);
}
struct double_comp_big {
double x;
int64_t i;
};
int compare_double_big( const void * l, const void * r )
{
const double * restrict _l= l;
const double * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct double_comp_big), compare_double_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_double_noidx_big(double* A, int64_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(double), compare_double_big);
}
struct float_comp {
float x;
int32_t i;
};
int compare_float( const void * l, const void * r )
{
const float * restrict _l= l;
const float * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct float_comp), compare_float);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_float_noidx(float* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(float), compare_float);
}
struct float_comp_big {
float x;
int64_t i;
};
int compare_float_big( const void * l, const void * r )
{
const float * restrict _l= l;
const float * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct float_comp_big), compare_float_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_float_noidx_big(float* A, int64_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(float), compare_float_big);
}
/* Generated C file:1 ends here */

169
src/utils/qsort.org Normal file
View File

@ -0,0 +1,169 @@
#+TITLE: Quick sort binding for Fortran
* C template
#+NAME: c_template
#+BEGIN_SRC c
struct TYPE_comp_big {
TYPE x;
int32_t i;
};
int compare_TYPE_big( const void * l, const void * r )
{
const TYPE * restrict _l= l;
const TYPE * restrict _r= r;
if( *_l > *_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<isize ; ++i) {
A[i].x = A_in[i];
A[i].i = iorder[i];
}
qsort( (void*) A, (size_t) isize, sizeof(struct TYPE_comp_big), compare_TYPE_big);
for (int i=0 ; i<isize ; ++i) {
A_in[i] = A[i].x;
iorder[i] = A[i].i;
}
free(A);
}
void qsort_TYPE_noidx_big(TYPE* A, int32_t isize) {
qsort( (void*) A, (size_t) isize, sizeof(TYPE), compare_TYPE_big);
}
#+END_SRC
* Fortran template
#+NAME:f_template
#+BEGIN_SRC f90
subroutine Lsort_big_c(A, iorder, isize) bind(C, name="qsort_TYPE_big")
use iso_c_binding
integer(c_int32_t), value :: isize
integer(c_int32_t) :: iorder(isize)
real (c_TYPE) :: A(isize)
end subroutine Lsort_big_c
subroutine Lsort_noidx_big_c(A, isize) bind(C, name="qsort_TYPE_noidx_big")
use iso_c_binding
integer(c_int32_t), value :: isize
real (c_TYPE) :: A(isize)
end subroutine Lsort_noidx_big_c
#+END_SRC
#+NAME:f_template2
#+BEGIN_SRC f90
subroutine Lsort_big(A, iorder, isize)
use qsort_module
use iso_c_binding
integer(c_int32_t) :: isize
integer(c_int32_t) :: iorder(isize)
real (c_TYPE) :: A(isize)
call Lsort_big_c(A, iorder, isize)
end subroutine Lsort_big
subroutine Lsort_noidx_big(A, isize)
use iso_c_binding
use qsort_module
integer(c_int32_t) :: isize
real (c_TYPE) :: A(isize)
call Lsort_noidx_big_c(A, isize)
end subroutine Lsort_noidx_big
#+END_SRC
* Python scripts for type replacements
#+NAME: replaced
#+begin_src python :results output :noweb yes
data = """
<<c_template>>
"""
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 = """
<<f_template>>
"""
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 = """
<<f_template2>>
"""
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 <stdlib.h>
#include <stdint.h>
<<replaced()>>
#+END_SRC
* Generated Fortran file
#+BEGIN_SRC f90 :tangle qsort_module.f90 :noweb yes
module qsort_module
use iso_c_binding
interface
<<replaced_f()>>
end interface
end module qsort_module
<<replaced_f2()>>
#+END_SRC

347
src/utils/qsort_module.f90 Normal file
View File

@ -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

2
src/utils/set_multiple_levels_omp.irp.f
View File

@ -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.)

695
src/utils/sort.irp.f
View File

@ -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 (j<k)
if ( x(j) < x(j+1) ) then
j=j+1
endif
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = shiftl(j,1)
else
j = k+1
endif
enddo
if (j==k) then
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = shiftl(j,1)
else
j = k+1
endif
endif
x(i) = xtemp
iorder(i) = i0
enddo
end subroutine heap_$Xsort
subroutine heap_$Xsort_big(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.
! This is a version for very large arrays where the indices need
! to be in integer*8 format
END_DOC
integer*8,intent(in) :: isize
$type,intent(inout) :: x(isize)
integer*8,intent(inout) :: iorder(isize)
integer*8 :: 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 (j<k)
if ( x(j) < x(j+1) ) then
j=j+1
endif
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = shiftl(j,1)
else
j = k+1
endif
enddo
if (j==k) then
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = shiftl(j,1)
else
j = k+1
endif
endif
x(i) = xtemp
iorder(i) = i0
enddo
end subroutine heap_$Xsort_big
subroutine sorted_$Xnumber(x,isize,n)
implicit none
@ -250,222 +32,6 @@ SUBST [ X, type ]
END_TEMPLATE
!---------------------- INTEL
IRP_IF INTEL
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize)
use intel
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
character, allocatable :: tmp(:)
if (isize < 2) return
call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
allocate(tmp(n))
call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp)
deallocate(tmp)
iorder(1:isize) = iorder(1:isize)+1
call $Xset_order(x,iorder,isize)
end
subroutine $Xsort_noidx(x,isize)
use intel
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
character, allocatable :: tmp(:)
if (isize < 2) return
call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
allocate(tmp(n))
call ippsSortRadixAscend_$ityp_I(x, isize, tmp)
deallocate(tmp)
end
SUBST [ X, type, ityp, n, ippsz ]
; real ; 32f ; 4 ; 13 ;;
i ; integer ; 32s ; 4 ; 11 ;;
i2 ; integer*2 ; 16s ; 2 ; 7 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize)
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
if (isize < 2) then
return
endif
! call sorted_$Xnumber(x,isize,n)
! if (isize == n) then
! return
! endif
if ( isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
! call heap_$Xsort(x,iorder,isize)
call quick_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
SUBST [ X, type ]
d ; double precision ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize)
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
if (isize < 2) then
return
endif
call sorted_$Xnumber(x,isize,n)
if (isize == n) then
return
endif
if ( isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
! call $Xradix_sort(x,iorder,isize,-1)
call quick_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
SUBST [ X, type ]
i8 ; integer*8 ;;
END_TEMPLATE
!---------------------- END INTEL
IRP_ELSE
!---------------------- NON-INTEL
BEGIN_TEMPLATE
subroutine $Xsort_noidx(x,isize)
implicit none
BEGIN_DOC
! Sort array x(isize).
END_DOC
integer,intent(in) :: isize
$type,intent(inout) :: x(isize)
integer, allocatable :: iorder(:)
integer :: i
allocate(iorder(isize))
do i=1,isize
iorder(i)=i
enddo
call $Xsort(x,iorder,isize)
deallocate(iorder)
end subroutine $Xsort_noidx
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8 ; integer*8 ;;
i2 ; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize)
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
if (isize < 2) then
return
endif
! call sorted_$Xnumber(x,isize,n)
! if (isize == n) then
! return
! endif
if ( isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
! call heap_$Xsort(x,iorder,isize)
call quick_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
SUBST [ X, type ]
; real ;;
d ; double precision ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize)
implicit none
BEGIN_DOC
! Sort array x(isize).
! 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 :: n
if (isize < 2) then
return
endif
call sorted_$Xnumber(x,isize,n)
if (isize == n) then
return
endif
if ( isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
! call $Xradix_sort(x,iorder,isize,-1)
call quick_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
SUBST [ X, type ]
i ; integer ;;
i8 ; integer*8 ;;
i2 ; integer*2 ;;
END_TEMPLATE
IRP_ENDIF
!---------------------- END NON-INTEL
BEGIN_TEMPLATE
subroutine $Xset_order(x,iorder,isize)
@ -491,47 +57,6 @@ BEGIN_TEMPLATE
deallocate(xtmp)
end
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8; integer*8 ;;
i2; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine insertion_$Xsort_big (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.
! This is a version for very large arrays where the indices need
! to be in integer*8 format
END_DOC
integer*8,intent(in) :: isize
$type,intent(inout) :: x(isize)
integer*8,intent(inout) :: iorder(isize)
$type :: xtmp
integer*8 :: i, i0, j, jmax
do i=2_8,isize
xtmp = x(i)
i0 = iorder(i)
j = i-1_8
do while (j>0_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

22
src/utils/units.irp.f Normal file
View File

@ -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

4
src/utils/util.irp.f
View File

@ -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

2
tests/bats/common.bats.sh
View File

@ -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

16
travis/compilation.sh
View File

@ -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

10
travis/configuration.sh
View File

@ -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

16
travis/testing.sh
View File

@ -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