LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Merge pull request #63 from QuantumPackage/bugfix

Browse Source 
Bugfix
This commit is contained in:
Anthony Scemama 2019-09-02 17:54:23 +02:00 committed by GitHub
commit 618d4ca8b3
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
41 changed files with 1359 additions and 424 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
INSTALL.rst
View File

@ -143,7 +143,7 @@ IRPF90
to Parameters (IRP) method.
* Download the latest version of IRPF90
here : `<https://github.com/scemama/irpf90/releases/latest>`_ and move
here : `<https://gitlab.com/scemama/irpf90/-/archive/v1.7.2/irpf90-v1.7.2.tar.gz>`_ and move
the downloaded archive in the :file:`${QP_ROOT}/external` directory
* Extract the archive and go into the :file:`irpf90-*` directory to run

1
REPLACE
View File

@ -834,3 +834,4 @@ qp_name potential_sr_c_alpha_ao_pbe --rename=potential_c_alpha_ao_sr_pbe
qp_name potential_sr_c_beta_ao_pbe --rename=potential_c_beta_ao_sr_pbe
qp_name potential_sr_xc_alpha_ao_pbe --rename=potential_xc_alpha_ao_sr_pbe
qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe
qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion

71
bin/qp_tunnel
View File

@ -1,71 +0,0 @@
#!/usr/bin/env python2
"""
Creates an ssh tunnel for using slaves on another network.
Launch a server on the front-end node of the cluster on which the master
process runs. Then start a client ont the front-end node of the distant
cluster.
Usage:
qp_tunnel server EZFIO_DIR
qp_tunnel client <address> EZFIO_DIR
Options:
-h --help
"""
import os
import sys
import zmq
try:
import qp_path
except ImportError:
print "source .quantum_package.rc"
raise
from docopt import docopt
from ezfio import ezfio
def get_address(filename):
with open(os.path.join(filename,'work','qp_run_address'),'r') as f:
a = f.readlines()[0].strip()
return a
def set_address(filename,address):
with open(os.path.join(filename,'work','qp_run_address'),'r') as f:
backup = f.readlines()
with open(os.path.join(filename,'work','qp_run_address'),'w') as f:
f.write('\n'.join([address]+backup))
def main_server(arguments,filename):
destination = get_address(filename)
print destination
def main_client(arguments,filename):
destination = arguments["<address>"]
print destination
def main(arguments):
"""Main function"""
print arguments
filename = arguments["EZFIO_DIR"]
if arguments["server"]:
return main_server(arguments, filename)
if arguments["client"]:
return main_client(arguments, filename)
if __name__ == '__main__':
ARGUMENTS = docopt(__doc__)
main(ARGUMENTS)

12
configure vendored
View File

@ -438,18 +438,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
fail
fi
BWRAP=$(find_exe bwrap)
if [[ ${BWRAP} = $(not_found) ]] ; then
error "Bubblewrap (bwrap) is not installed."
fail
fi
LIBCAP=$(find_lib -lcap)
if [[ ${LIBCAP} = $(not_found) ]] ; then
error "Libcap (libcap) is not installed."
fail
fi
BWRAP=$(find_exe bwrap)
if [[ ${BWRAP} = $(not_found) ]] ; then
error "Bubblewrap (bwrap) is not installed."
fail
fi
OPAM=$(find_exe opam)
if [[ ${OPAM} = $(not_found) ]] ; then
error "OPAM (ocaml) package manager is not installed."

BIN
data/qp2.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 5.9 MiB

After

Width:  |  Height:  |  Size: 351 KiB

BIN
data/qp2_hd.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 5.9 MiB

89
docs/source/users_guide/qp_tunnel.rst Normal file
View File

@ -0,0 +1,89 @@
.. _qp_tunnel:
=========
qp_tunnel
=========
.. TODO
.. program:: qp_tunnel
Establishes a tunnel to allow communications between machines within
different networks, for example multiple MPI slave jobs running on
different clusters.
Usage
-----
.. code:: bash
qp_tunnel [-g] (ADDRESS|EZFIO_DIR)
``EZFIO_DIR`` is the name of the |EZFIO| directory containing the data,
and ``ADDRESS`` is the address of another tunnel.
.. option:: -h, --help
Displays the help message
.. option:: -g, --get-input
Download the EZFIO directory from the remote instance of qp_tunnel.
Example
-------
.. code:: text
+-------------------+ +------------------+
| | | |
| N1_1 N1_2 N1_3 | | N2_1 N2_2 N2_3 |
| | | | | | | | | |
| +----+----+ | | +----+----+ |
| | | | | |
| C1 F1 | | F2 C2 |
| +---------=----=--------+ |
| | | |
+-------------------+ +------------------+
Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH
on a front-end named respectively F1 and F2. Groups of nodes N1 and N2 have
been reserved by the batch scheduling system on both clusters.
Each node in N1 is on the same network as the other nodes of N1, but they
can't access the network on which the nodes of N2 are.
1) Start a parallel simulation on the cluster C1, running on nodes N1.
We assume that there is a shared file system, such that F1 can access
the EZFIO directory. We also assume that F1 can communicate with the
nodes of N1.
2) Run a tunnel on the front-end F1 and keep it running:
.. code:: bash
me@f1 $ qp_tunnel my_directory.ezfio
Connect to:
tcp://31.122.230.47:42379
Ready
3) On the front-end F2, run another instance connecting to the other one,
which will fetch the |EZFIO| directory:
.. code:: bash
me@f2 $ qp_tunnel --get-input tcp://31.122.230.47:42379
Connect to:
tcp://31.122.209.139:42379
Communication [ OK ]
Getting input... my_directory.ezfio ...done
Ready
4) Keep the tunnel running, and you can now run a slave simulation within the
nodes N2.

139
man/qp_tunnel.1 Normal file
View File

@ -0,0 +1,139 @@
.\" Man page generated from reStructuredText.
.
.TH "QP_TUNNEL" "1" "Jun 15, 2019" "2.0" "Quantum Package"
.SH NAME
qp_tunnel \- | Quantum Package >
.
.nr rst2man-indent-level 0
.
.de1 rstReportMargin
\\$1 \\n[an-margin]
level \\n[rst2man-indent-level]
level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
-
\\n[rst2man-indent0]
\\n[rst2man-indent1]
\\n[rst2man-indent2]
..
.de1 INDENT
.\" .rstReportMargin pre:
. RS \\$1
. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
. nr rst2man-indent-level +1
.\" .rstReportMargin post:
..
.de UNINDENT
. RE
.\" indent \\n[an-margin]
.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
.nr rst2man-indent-level -1
.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
..
.sp
Establishes a tunnel to allow communications between machines within
different networks, for example multiple MPI slave jobs running on
different clusters.
.SH USAGE
.INDENT 0.0
.INDENT 3.5
.sp
.nf
.ft C
qp_tunnel [\-g] (ADDRESS|EZFIO_DIR)
.ft P
.fi
.UNINDENT
.UNINDENT
.sp
\fBEZFIO_DIR\fP is the name of the \fI\%EZFIO\fP directory containing the data,
and \fBADDRESS\fP is the address of another tunnel.
.INDENT 0.0
.TP
.B \-h, \-\-help
Displays the help message
.UNINDENT
.INDENT 0.0
.TP
.B \-g, \-\-get\-input
Download the EZFIO directory from the remote instance of qp_tunnel.
.UNINDENT
.SH EXAMPLE
.INDENT 0.0
.INDENT 3.5
.sp
.nf
.ft C
+\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+
| | | |
| N1_1 N1_2 N1_3 | | N2_1 N2_2 N2_3 |
| | | | | | | | | |
| +\-\-\-\-+\-\-\-\-+ | | +\-\-\-\-+\-\-\-\-+ |
| | | | | |
| C1 F1 | | F2 C2 |
| +\-\-\-\-\-\-\-\-\-=\-\-\-\-=\-\-\-\-\-\-\-\-+ |
| | | |
+\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+
.ft P
.fi
.UNINDENT
.UNINDENT
.sp
Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH
on a front\-end named respectively F1 and F2. Groups of nodes N1 and N2 have
been reserved by the batch scheduling system on both clusters.
Each node in N1 is on the same network as the other nodes of N1, but they
can’t access the network on which the nodes of N2 are.
.INDENT 0.0
.IP 1. 3
Start a parallel simulation on the cluster C1, running on nodes N1.
We assume that there is a shared file system, such that F1 can access
the EZFIO directory. We also assume that F1 can communicate with the
nodes of N1.
.IP 2. 3
Run a tunnel on the front\-end F1 and keep it running:
.UNINDENT
.INDENT 0.0
.INDENT 3.5
.sp
.nf
.ft C
me@f1 $ qp_tunnel my_directory.ezfio
Connect to:
tcp://31.122.230.47:42379
Ready
.ft P
.fi
.UNINDENT
.UNINDENT
.INDENT 0.0
.IP 3. 3
On the front\-end F2, run another instance connecting to the other one,
which will fetch the \fI\%EZFIO\fP directory:
.UNINDENT
.INDENT 0.0
.INDENT 3.5
.sp
.nf
.ft C
me@f2 $ qp_tunnel \-\-get\-input tcp://31.122.230.47:42379
Connect to:
tcp://31.122.209.139:42379
Communication [ OK ]
Getting input... my_directory.ezfio ...done
Ready
.ft P
.fi
.UNINDENT
.UNINDENT
.INDENT 0.0
.IP 4. 3
Keep the tunnel running, and you can now run a slave simulation within the
nodes N2.
.UNINDENT
.SH AUTHOR
A. Scemama, E. Giner
.SH COPYRIGHT
2019, A. Scemama, E. Giner
.\" Generated by docutils manpage writer.
.

2
man/rotate_mos.1
View File

@ -34,7 +34,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
.INDENT 3.5
Rotates molecular orbitals i and j by combining them as
$1/sqrt{2} ( phi_i + phi_j )$ and
$1/sqrt{2} ( phi_i \- phi_j )$.
$1/sqrt{2} ( phi_i - phi_j )$.
.sp
Needs:
.INDENT 0.0

55
man/test.1 Normal file
View File

@ -0,0 +1,55 @@
.\" Man page generated from reStructuredText.
.
.TH "TEST" "1" "Jun 15, 2019" "2.0" "Quantum Package"
.SH NAME
test \- | Quantum Package >
.
.nr rst2man-indent-level 0
.
.de1 rstReportMargin
\\$1 \\n[an-margin]
level \\n[rst2man-indent-level]
level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
-
\\n[rst2man-indent0]
\\n[rst2man-indent1]
\\n[rst2man-indent2]
..
.de1 INDENT
.\" .rstReportMargin pre:
. RS \\$1
. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
. nr rst2man-indent-level +1
.\" .rstReportMargin post:
..
.de UNINDENT
. RE
.\" indent \\n[an-margin]
.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
.nr rst2man-indent-level -1
.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
..
.INDENT 0.0
.INDENT 3.5
Calls:
.INDENT 0.0
.INDENT 2.0
.IP \(bu 2
\fBtwo_e_integrals_index()\fP
.UNINDENT
.INDENT 2.0
.IP \(bu 2
\fBtwo_e_integrals_index_reverse()\fP
.UNINDENT
.INDENT 2.0
.UNINDENT
.UNINDENT
.UNINDENT
.UNINDENT
.SH AUTHOR
A. Scemama, E. Giner
.SH COPYRIGHT
2019, A. Scemama, E. Giner
.\" Generated by docutils manpage writer.
.

2
ocaml/.merlin
View File

@ -1,4 +1,4 @@
PKG core ZMQ cryptokit
PKG core zmq cryptokit
B _build/

214
ocaml/Bitlist.ml
View File

@ -7,82 +7,61 @@ Type for bits strings
list of Bits
*)
type t = Bit.t list
type t = int64 array
let n_int = Array.length
(* Create a zero bit list *)
let zero n_int =
Array.make (N_int_number.to_int n_int) 0L
(* String representation *)
let to_string b =
let rec do_work accu = function
| [] -> accu
| head :: tail ->
let new_accu = (Bit.to_string head) ^ accu
in do_work new_accu tail
let int64_to_string x =
String.init 64 (fun i ->
if Int64.logand x @@ Int64.shift_left 1L i <> 0L then
'+'
else
'-')
in
do_work "" b
Array.map int64_to_string b
|> Array.to_list
|> String.concat ""
let of_string ?(zero='0') ?(one='1') s =
List.init (String.length s) (String.get s)
|> List.rev_map ( fun c ->
if (c = zero) then Bit.Zero
else if (c = one) then Bit.One
else (failwith ("Error in bitstring ") ) )
let n_int = ( (String.length s - 1) lsr 6 ) + 1 in
let result = Array.make n_int 0L in
String.iteri (fun i c ->
if c = one then
begin
let iint = i lsr 6 in (* i / 64 *)
let k = i - (iint lsl 6) in
result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
end) s;
result
let of_string_mp s =
List.init (String.length s) (String.get s)
|> List.rev_map (function
| '-' -> Bit.Zero
| '+' -> Bit.One
| _ -> failwith ("Error in bitstring ") )
let of_string_mp = of_string ~zero:'-' ~one:'+'
(* Create a bit list from an int64 *)
let of_int64 i =
let rec do_work accu = function
| 0L -> Bit.Zero :: accu |> List.rev
| 1L -> Bit.One :: accu |> List.rev
| i ->
let b =
match (Int64.logand i 1L ) with
| 0L -> Bit.Zero
| 1L -> Bit.One
| _ -> raise (Failure "i land 1 not in (0,1)")
in
do_work (b :: accu) (Int64.shift_right_logical i 1)
in
let adjust_length result =
let rec do_work accu = function
| 64 -> List.rev accu
| i when i>64 -> raise (Failure "Error in of_int64 > 64")
| i when i<0 -> raise (Failure "Error in of_int64 < 0")
| i -> do_work (Bit.Zero :: accu) (i+1)
in
do_work (List.rev result) (List.length result)
in
adjust_length (do_work [] i)
let of_int64 i = [| i |]
(* Create an int64 from a bit list *)
let to_int64 l =
assert ( (List.length l) <= 64) ;
let rec do_work accu = function
| [] -> accu
| Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail
| Bit.One::tail -> do_work Int64.(logor one (shift_left accu 1)) tail
in do_work Int64.zero (List.rev l)
let to_int64 = function
| [| i |] -> i
| _ -> failwith "N_int > 1"
(* Create a bit list from a list of int64 *)
let of_int64_list l =
List.map of_int64 l
|> List.concat
(* Create a bit list from an array of int64 *)
let of_int64_array l =
Array.map of_int64 l
|> Array.to_list
|> List.concat
external of_int64_array : int64 array -> t = "%identity"
external to_int64_array : t -> int64 array = "%identity"
(* Create a bit list from a list of int64 *)
let of_int64_list l =
Array.of_list l |> of_int64_array
(* Compute n_int *)
@ -91,101 +70,64 @@ let n_int_of_mo_num mo_num =
N_int_number.of_int ( (mo_num-1)/bit_kind_size + 1 )
(* Create a zero bit list *)
let zero n_int =
let n_int = N_int_number.to_int n_int in
let a = Array.init n_int (fun i-> 0L) in
of_int64_list ( Array.to_list a )
(* Create an int64 list from a bit list *)
let to_int64_list l =
let rec do_work accu buf counter = function
| [] ->
begin
match buf with
| [] -> accu
| _ -> (List.rev buf)::accu
end
| i::tail ->
if (counter < 64) then
do_work accu (i::buf) (counter+1) tail
else
do_work ( (List.rev (i::buf))::accu) [] 1 tail
in
let l = do_work [] [] 1 l
in
List.rev_map to_int64 l
to_int64_array l |> Array.to_list
(* Create an array of int64 from a bit list *)
let to_int64_array l =
to_int64_list l
|> Array.of_list
(* Create a bit list from a list of MO indices *)
let of_mo_number_list n_int l =
let n_int = N_int_number.to_int n_int in
let length = n_int*64 in
let a = Array.make length (Bit.Zero) in
List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
Array.to_list a
let result = zero n_int in
List.iter (fun j ->
let i = (MO_number.to_int j) - 1 in
let iint = i lsr 6 in (* i / 64 *)
let k = i - (iint lsl 6) in
result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
) l;
result
let to_mo_number_list l =
let a = Array.of_list l in
let mo_num = MO_number.get_max () in
let rec do_work accu = function
| 0 -> accu
| i ->
begin
let new_accu =
match a.(i-1) with
| Bit.One -> (MO_number.of_int ~max:mo_num i)::accu
| Bit.Zero -> accu
in
do_work new_accu (i-1)
end
let rec aux_one x shift accu = function
| -1 -> accu
| i -> if Int64.logand x (Int64.shift_left 1L i) <> 0L then
aux_one x shift ( (i+shift) ::accu) (i-1)
else
aux_one x shift accu (i-1)
in
do_work [] (List.length l)
Array.mapi (fun i x ->
let shift = (i lsr 6) lsl 6 + 1 in
aux_one x shift [] 63
) l
|> Array.to_list
|> List.concat
|> List.map MO_number.of_int
(* logical operations on bit_list *)
let logical_operator2 op a b =
let rec do_work_binary result a b =
match a, b with
| [], [] -> result
| [], _ | _ , [] -> raise (Failure "Lists should have same length")
| (ha::ta), (hb::tb) ->
let newbit = op ha hb
in do_work_binary (newbit::result) ta tb
let and_operator a b = Array.map2 Int64.logand a b
let xor_operator a b = Array.map2 Int64.logxor a b
let or_operator a b = Array.map2 Int64.logor a b
let not_operator b = Array.map Int64.lognot b
let pop_sign =
let mask =
(Int64.pred (Int64.shift_left 1L 63))
in
List.rev (do_work_binary [] a b)
let logical_operator1 op b =
let rec do_work_unary result b =
match b with
| [] -> result
| (hb::tb) ->
let newbit = op hb
in do_work_unary (newbit::result) tb
in
List.rev (do_work_unary [] b)
let and_operator a b = logical_operator2 Bit.and_operator a b
let xor_operator a b = logical_operator2 Bit.xor_operator a b
let or_operator a b = logical_operator2 Bit.or_operator a b
let not_operator b = logical_operator1 Bit.not_operator b
fun x -> Int64.logand mask x
let popcnt b =
List.fold_left (fun accu -> function
| Bit.One -> accu+1
| Bit.Zero -> accu
) 0 b
Array.fold_left (fun accu x ->
if x >= 0L then
accu + (Z.popcount @@ Z.of_int64 x)
else
accu + 1 + (Z.popcount @@ Z.of_int64 (pop_sign x))
) 0 b

2
ocaml/Bitlist.mli
View File

@ -1,4 +1,4 @@
type t = Bit.t list
type t
(** The zero bit list *)
val zero : Qptypes.N_int_number.t -> t

31
ocaml/Determinant.ml
View File

@ -25,19 +25,6 @@ let to_bitlist_couple x =
in (xa,xb)
let bitlist_to_string ~mo_num x =
let len =
MO_number.to_int mo_num
in
let s =
List.map (function
| Bit.Zero -> "-"
| Bit.One -> "+"
) x
|> String.concat ""
in
String.sub s 0 len
let of_int64_array ~n_int ~alpha ~beta x =
@ -48,37 +35,29 @@ let of_int64_array ~n_int ~alpha ~beta x =
in
if ( (Bitlist.popcnt a) <> alpha) then
begin
let mo_num = MO_number.get_max () in
let mo_num = MO_number.of_int mo_num ~max:mo_num in
failwith (Printf.sprintf "Expected %d electrons in alpha determinant
%s" alpha (bitlist_to_string ~mo_num:mo_num a) )
%s" alpha (Bitlist.to_string a) )
end;
if ( (Bitlist.popcnt b) <> beta ) then
begin
let mo_num = MO_number.get_max () in
let mo_num = MO_number.of_int mo_num ~max:mo_num in
failwith (Printf.sprintf "Expected %d electrons in beta determinant
%s" beta (bitlist_to_string ~mo_num:mo_num b) )
%s" beta (Bitlist.to_string b) )
end;
x
let of_bitlist_couple ?n_int ~alpha ~beta (xa,xb) =
let of_bitlist_couple ~n_int ~alpha ~beta (xa,xb) =
let ba, bb =
Bitlist.to_int64_array xa ,
Bitlist.to_int64_array xb
and n_int =
match n_int with
| Some x -> x
| None -> Bitlist.n_int_of_mo_num (List.length xa)
in
of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb])
let to_string ~mo_num x =
let (xa,xb) = to_bitlist_couple x in
[ " " ; bitlist_to_string ~mo_num xa ; "\n" ;
" " ; bitlist_to_string ~mo_num xb ]
[ " " ; Bitlist.to_string xa ; "\n" ;
" " ; Bitlist.to_string xb ]
|> String.concat ""

2
ocaml/Determinant.mli
View File

@ -24,7 +24,7 @@ val to_alpha_beta : t -> (int64 array)*(int64 array)
val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
(** Create from a bit list *)
val of_bitlist_couple : ?n_int:Qptypes.N_int_number.t ->
val of_bitlist_couple : n_int:Qptypes.N_int_number.t ->
alpha:Qptypes.Elec_alpha_number.t ->
beta:Qptypes.Elec_beta_number.t ->
Bitlist.t * Bitlist.t -> t

133
ocaml/Input_determinants_by_hand.ml
View File

@ -7,14 +7,14 @@ module Determinants_by_hand : sig
{ n_int : N_int_number.t;
bit_kind : Bit_kind.t;
n_det : Det_number.t;
n_det_qp_edit : Det_number.t;
n_states : States_number.t;
expected_s2 : Positive_float.t;
psi_coef : Det_coef.t array;
psi_det : Determinant.t array;
state_average_weight : Positive_float.t array;
} [@@deriving sexp]
val read : unit -> t
val read_maybe : unit -> t option
val read : ?full:bool -> unit -> t option
val write : t -> unit
val to_string : t -> string
val to_rst : t -> Rst_string.t
@ -28,6 +28,7 @@ end = struct
{ n_int : N_int_number.t;
bit_kind : Bit_kind.t;
n_det : Det_number.t;
n_det_qp_edit : Det_number.t;
n_states : States_number.t;
expected_s2 : Positive_float.t;
psi_coef : Det_coef.t array;
@ -38,8 +39,6 @@ end = struct
let get_default = Qpackage.get_ezfio_default "determinants";;
let n_det_read_max = 10_000 ;;
let read_n_int () =
if not (Ezfio.has_determinants_n_int()) then
Ezfio.get_mo_basis_mo_num ()
@ -80,11 +79,27 @@ end = struct
|> Det_number.of_int
;;
let read_n_det_qp_edit () =
if not (Ezfio.has_determinants_n_det_qp_edit ()) then
begin
let n_det = read_n_det () |> Det_number.to_int in
Ezfio.set_determinants_n_det_qp_edit n_det
end;
Ezfio.get_determinants_n_det_qp_edit ()
|> Det_number.of_int
;;
let write_n_det n =
Det_number.to_int n
|> Ezfio.set_determinants_n_det
;;
let write_n_det_qp_edit n =
let n_det = read_n_det () |> Det_number.to_int in
min n_det (Det_number.to_int n)
|> Ezfio.set_determinants_n_det_qp_edit
;;
let read_n_states () =
if not (Ezfio.has_determinants_n_states ()) then
Ezfio.set_determinants_n_states 1
@ -178,7 +193,7 @@ end = struct
|> Ezfio.set_determinants_expected_s2
;;
let read_psi_coef () =
let read_psi_coef ~read_only () =
if not (Ezfio.has_determinants_psi_coef ()) then
begin
let n_states =
@ -189,7 +204,12 @@ end = struct
~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. ))
|> Ezfio.set_determinants_psi_coef
end;
Ezfio.get_determinants_psi_coef ()
begin
if read_only then
Ezfio.get_determinants_psi_coef_qp_edit ()
else
Ezfio.get_determinants_psi_coef ()
end
|> Ezfio.flattened_ezfio
|> Array.map Det_coef.of_float
;;
@ -202,12 +222,15 @@ end = struct
and n_states =
States_number.to_int n_states
in
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
|> Ezfio.set_determinants_psi_coef
let r =
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
in
Ezfio.set_determinants_psi_coef r;
Ezfio.set_determinants_psi_coef_qp_edit r
;;
let read_psi_det () =
let read_psi_det ~read_only () =
let n_int = read_n_int ()
and n_alpha = Ezfio.get_electrons_elec_alpha_num ()
|> Elec_alpha_number.of_int
@ -233,13 +256,21 @@ end = struct
|> Ezfio.set_determinants_psi_det ;
end ;
let n_int = N_int_number.to_int n_int in
let psi_det_array = Ezfio.get_determinants_psi_det () in
let psi_det_array =
if read_only then
Ezfio.get_determinants_psi_det_qp_edit ()
else
Ezfio.get_determinants_psi_det ()
in
let dim = psi_det_array.Ezfio.dim
and data = Ezfio.flattened_ezfio psi_det_array
in
assert (n_int = dim.(0));
assert (dim.(1) = 2);
assert (dim.(2) = (Det_number.to_int (read_n_det ())));
if read_only then
assert (dim.(2) = (Det_number.to_int (read_n_det_qp_edit ())))
else
assert (dim.(2) = (Det_number.to_int (read_n_det ())));
Array.init dim.(2) (fun i ->
Array.sub data (2*n_int*i) (2*n_int) )
|> Array.map (Determinant.of_int64_array
@ -252,40 +283,45 @@ end = struct
|> Array.concat
|> Array.to_list
in
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
|> Ezfio.set_determinants_psi_det
let r =
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
in
Ezfio.set_determinants_psi_det r;
Ezfio.set_determinants_psi_det_qp_edit r
;;
let read () =
let read ?(full=true) () =
let n_det_qp_edit = read_n_det_qp_edit () in
let n_det = read_n_det () in
let read_only =
if full then false else n_det_qp_edit <> n_det
in
if (Ezfio.has_mo_basis_mo_num ()) then
try
Some
{ n_int = read_n_int () ;
bit_kind = read_bit_kind () ;
n_det = read_n_det () ;
n_det_qp_edit = read_n_det_qp_edit () ;
expected_s2 = read_expected_s2 () ;
psi_coef = read_psi_coef () ;
psi_det = read_psi_det () ;
psi_coef = read_psi_coef ~read_only () ;
psi_det = read_psi_det ~read_only () ;
n_states = read_n_states () ;
state_average_weight = read_state_average_weight () ;
}
with _ -> None
else
failwith "No molecular orbitals, so no determinants"
;;
let read_maybe () =
let n_det =
read_n_det ()
in
if ( (Det_number.to_int n_det) < n_det_read_max ) then
try Some (read ()) with
| Failure _ -> None
else
(* No molecular orbitals, so no determinants *)
None
;;
let write { n_int ;
bit_kind ;
n_det ;
n_det_qp_edit ;
expected_s2 ;
psi_coef ;
psi_det ;
@ -297,9 +333,13 @@ end = struct
write_n_det n_det;
write_n_states n_states;
write_expected_s2 expected_s2;
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
write_state_average_weight state_average_weight;
if n_det <= n_det_qp_edit then
begin
write_n_det_qp_edit n_det;
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
write_psi_det ~n_int:n_int ~n_det:n_det psi_det
end;
write_state_average_weight state_average_weight
;;
@ -316,11 +356,13 @@ end = struct
|> States_number.to_int
and ndet =
Det_number.to_int b.n_det
and ndet_qp_edit =
Det_number.to_int b.n_det_qp_edit
in
let coefs_string i =
Array.init nstates (fun j ->
let ishift =
j*ndet
j*ndet_qp_edit
in
if (ishift < Array.length b.psi_coef) then
b.psi_coef.(i+ishift)
@ -331,7 +373,7 @@ end = struct
)
|> Array.to_list |> String.concat "\t"
in
Array.init ndet (fun i ->
Array.init ndet_qp_edit (fun i ->
Printf.sprintf " %s\n%s\n"
(coefs_string i)
(Determinant.to_string ~mo_num:mo_num b.psi_det.(i)
@ -472,6 +514,7 @@ psi_det = %s
(* Handle determinants *)
let psi_det =
let n_int = N_int_number.of_int @@ (MO_number.get_max () - 1) / 64 + 1 in
let n_alpha = Ezfio.get_electrons_elec_alpha_num ()
|> Elec_alpha_number.of_int
and n_beta = Ezfio.get_electrons_elec_beta_num ()
@ -483,8 +526,8 @@ psi_det = %s
begin
let newdet =
(Bitlist.of_string ~zero:'-' ~one:'+' alpha ,
Bitlist.of_string ~zero:'-' ~one:'+' beta)
|> Determinant.of_bitlist_couple ~alpha:n_alpha ~beta:n_beta
Bitlist.of_string ~zero:'-' ~one:'+' beta)
|> Determinant.of_bitlist_couple ~n_int ~alpha:n_alpha ~beta:n_beta
|> Determinant.sexp_of_t
|> Sexplib.Sexp.to_string
in
@ -492,9 +535,6 @@ psi_det = %s
end
| _::tail -> read_dets accu tail
in
let dets =
List.map String_ext.rev dets
in
let a =
read_dets [] dets
|> String.concat ""
@ -510,9 +550,11 @@ psi_det = %s
Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
and n_states =
Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ())
and n_det_qp_edit =
Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ())
in
let s =
String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det]
String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ]
in
@ -527,7 +569,9 @@ psi_det = %s
Det_number.to_int n_det_new
in
let det =
read ()
match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in
let n_det_old, n_states =
Det_number.to_int det.n_det,
@ -558,7 +602,9 @@ psi_det = %s
let extract_state istate =
Printf.printf "Extracting state %d\n" (States_number.to_int istate);
let det =
read ()
match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in
let n_det, n_states =
Det_number.to_int det.n_det,
@ -588,7 +634,9 @@ psi_det = %s
let extract_states range =
Printf.printf "Extracting states %s\n" (Range.to_string range);
let det =
read ()
match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in
let n_det, n_states =
Det_number.to_int det.n_det,
@ -614,7 +662,8 @@ psi_det = %s
j*n_det
in
for i=0 to (n_det-1) do
det.psi_coef.(!state_shift+i) <- det.psi_coef.(i+ishift)
det.psi_coef.(!state_shift+i) <-
det.psi_coef.(i+ishift)
done
end;
state_shift := !state_shift + n_det

11
ocaml/Input_mo_basis.ml
View File

@ -65,8 +65,15 @@ end = struct
let read_mo_num () =
Ezfio.get_mo_basis_mo_num ()
|> MO_number.of_int
let elec_alpha_num =
Ezfio.get_electrons_elec_alpha_num ()
in
let result =
Ezfio.get_mo_basis_mo_num ()
in
if result < elec_alpha_num then
failwith "More alpha electrons than MOs";
MO_number.of_int result
let read_mo_class () =

2
ocaml/Makefile
View File

@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m
$(QP_ROOT)/ocaml/element_create_db.byte
external_libs:
opam install cryptokit core
opam install cryptokit sexplib
qpackage.odocl: $(MLIFILES)
ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl

2
ocaml/_tags
View File

@ -1,4 +1,4 @@
true: package(cryptokit,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
true: thread
false: profile
<*byte> : linkdep(c_bindings.o), custom

9
ocaml/qp_print_basis.ml
View File

@ -44,9 +44,12 @@ let psi_det () =
let psi_det =
Input.Determinants_by_hand.read ()
in
Input.Determinants_by_hand.to_rst psi_det
|> Rst_string.to_string
|> print_endline
match psi_det with
| Some psi_det ->
Input.Determinants_by_hand.to_rst psi_det
|> Rst_string.to_string
|> print_endline
| None -> ()

11
ocaml/qp_set_mo_class.ml
View File

@ -112,9 +112,9 @@ let set ~core ~inact ~act ~virt ~del =
and av = Excitation.create_single act virt
in
let single_excitations = [ ia ; aa ; av ]
|> List.map (fun x ->
|> List.map (fun z ->
let open Excitation in
match x with
match z with
| Single (x,y) ->
( MO_class.to_bitlist n_int (Hole.to_mo_class x),
MO_class.to_bitlist n_int (Particle.to_mo_class y) )
@ -187,9 +187,10 @@ let set ~core ~inact ~act ~virt ~del =
match aa with
| Double _ -> assert false
| Single (x,y) ->
( MO_class.to_bitlist n_int (Hole.to_mo_class x) ) @
( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
|> Bitlist.to_int64_list
Bitlist.to_int64_list
( MO_class.to_bitlist n_int ( Hole.to_mo_class x) ) @
Bitlist.to_int64_list
( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
in
Ezfio.set_bitmasks_n_mask_cas 1;
Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result

470
ocaml/qp_tunnel.ml Normal file
View File

@ -0,0 +1,470 @@
open Qputils
open Qptypes
type ezfio_or_address = EZFIO of string | ADDRESS of string
type req_or_sub = REQ | SUB
let localport = 42379
let in_time_sum = ref 1.e-9
and in_size_sum = ref 0.
let () =
let open Command_line in
begin
"Creates an ssh tunnel for using slaves on another network. Launch a server on the front-end node of the cluster on which the master process runs. Then start a client ont the front-end node of the distant cluster."
|> set_footer_doc ;
[ { short='g' ; long="get-input" ; opt=Optional ;
doc="Downloads the EZFIO directory." ;
arg=Without_arg; } ;
{ short='v' ; long="verbose" ; opt=Optional ;
doc="Prints the transfer speed." ;
arg=Without_arg; } ;
anonymous
"(EZFIO_DIR|ADDRESS)"
Mandatory
"EZFIO directory or address.";
] |> set_specs
end;
let arg =
let x =
match Command_line.anon_args () with
| [x] -> x
| _ -> begin
Command_line.help () ;
failwith "EZFIO_FILE or ADDRESS is missing"
end
in
if Sys.file_exists x && Sys.is_directory x then
EZFIO x
else
ADDRESS x
in
let verbose =
Command_line.get_bool "verbose"
in
let localhost =
Lazy.force TaskServer.ip_address
in
let long_address =
match arg with
| ADDRESS x -> x
| EZFIO x ->
let ic =
Filename.concat (Qpackage.ezfio_work x) "qp_run_address"
|> open_in
in
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 @@
Printf.sprintf "%s : Malformed address" long_address
in
let zmq_context =
Zmq.Context.create ()
in
(** Check availability of the ports *)
let localport =
let dummy_socket =
Zmq.Socket.create zmq_context Zmq.Socket.rep
in
let rec try_new_port port_number =
try
List.iter (fun i ->
let address =
Printf.sprintf "tcp://%s:%d" localhost (port_number+i)
in
Zmq.Socket.bind dummy_socket address;
Zmq.Socket.unbind dummy_socket address
) [ 0;1;2;3;4;5;6;7;8;9 ] ;
port_number
with
| Unix.Unix_error _ -> try_new_port (port_number+100)
in
let result =
try_new_port localport
in
Zmq.Socket.close dummy_socket;
result
in
let create_socket sock_type bind_or_connect addr =
let socket =
Zmq.Socket.create zmq_context sock_type
in
let () =
try
bind_or_connect socket addr
with
| _ -> failwith @@
Printf.sprintf "Unable to establish connection to %s." addr
in
socket
in
(* Handle termination *)
let run_status = ref true in
let handler =
Sys.Signal_handle (fun signum ->
run_status := false;
Sys.set_signal signum Sys.Signal_default
)
in
Sys.set_signal Sys.sigusr1 handler;
Sys.set_signal Sys.sigint handler;
let new_thread req_or_sub 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 =
if verbose then
begin
match req_or_sub with
| REQ -> (fun () ->
let msg =
Zmq.Socket.recv_all socket_in
in
let t0 = Unix.gettimeofday () in
Zmq.Socket.send_all socket_out msg;
let in_size =
float_of_int ( List.fold_left (fun accu x -> accu + String.length x) 0 msg )
/. 8192. /. 1024.
in
let msg =
Zmq.Socket.recv_all socket_out
in
let t1 = Unix.gettimeofday () in
Zmq.Socket.send_all socket_in msg;
let in_time = t1 -. t0 in
in_time_sum := !in_time_sum +. in_time;
in_size_sum := !in_size_sum +. in_size;
Printf.printf " %16.2f MiB/s -- %16.2f MiB/s\n%!" (in_size /. in_time) (!in_size_sum /. !in_time_sum);
)
| SUB -> (fun () ->
Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out)
end
else
begin
match req_or_sub with
| REQ -> (fun () ->
Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out;
Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in )
| SUB -> (fun () ->
Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out)
end
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)
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 () -> () )
in
let pollitem =
Zmq.Poll.mask_of
[| (socket_in, Zmq.Poll.In) ; (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 ocaml_thread =
let addr_out =
Printf.sprintf "tcp:%s:%d" address port
in
let addr_in =
Printf.sprintf "tcp://*:%d" localport
in
let f () =
new_thread REQ addr_in addr_out
in
(Thread.create f) ()
in
Printf.printf "Connect to:\ntcp://%s:%d\n%!" localhost localport;
let fortran_thread =
let addr_out =
Printf.sprintf "tcp:%s:%d" address (port+2)
in
let addr_in =
Printf.sprintf "tcp://*:%d" (localport+2)
in
let f () =
new_thread REQ addr_in addr_out
in
(Thread.create f) ()
in
let pub_thread =
let addr_in =
Printf.sprintf "tcp:%s:%d" address (port+1)
in
let addr_out =
Printf.sprintf "tcp://*:%d" (localport+1)
in
let f () =
new_thread SUB addr_in addr_out
in
(Thread.create f) ()
in
let input_thread =
let f () =
let addr_out =
match arg with
| EZFIO _ -> None
| ADDRESS _ -> Some (
Printf.sprintf "tcp:%s:%d" address (port+9) )
in
let addr_in =
Printf.sprintf "tcp://*:%d" (localport+9)
in
let socket_in =
create_socket Zmq.Socket.rep Zmq.Socket.bind addr_in
in
let socket_out =
match addr_out with
| Some addr_out -> Some (
create_socket Zmq.Socket.req Zmq.Socket.connect addr_out)
| None -> None
in
let temp_file =
Filename.temp_file "qp_tunnel" ".tar.gz"
in
let get_ezfio_filename () =
match arg with
| EZFIO x -> x
| ADDRESS _ ->
begin
match socket_out with
| None -> assert false
| Some socket_out -> (
Zmq.Socket.send socket_out "get_ezfio_filename" ;
Zmq.Socket.recv socket_out
)
end
in
let get_input () =
match arg with
| EZFIO x ->
begin
Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x
|> Sys.command |> ignore;
let fd =
Unix.openfile temp_file [Unix.O_RDONLY] 0o640
in
let len =
Unix.lseek fd 0 Unix.SEEK_END
in
ignore @@ Unix.lseek fd 0 Unix.SEEK_SET ;
let bstr =
Unix.map_file fd Bigarray.char
Bigarray.c_layout false [| len |]
|> Bigarray.array1_of_genarray
in
let result =
String.init len (fun i -> bstr.{i}) ;
in
Unix.close fd;
Sys.remove temp_file;
result
end
| ADDRESS _ ->
begin
match socket_out with
| None -> assert false
| Some socket_out -> (
Zmq.Socket.send socket_out "get_input" ;
Zmq.Socket.recv socket_out
)
end
in
let () =
match socket_out with
| None -> ()
| Some socket_out ->
Zmq.Socket.send socket_out "test";
Printf.printf "Communication [ %s ]\n%!" (Zmq.Socket.recv socket_out);
in
(* Download input if asked *)
if Command_line.get_bool "get-input" then
begin
match arg with
| EZFIO _ -> ()
| ADDRESS _ ->
begin
Printf.printf "Getting input... %!";
let ezfio_filename =
get_ezfio_filename ()
in
Printf.printf "%s%!" ezfio_filename;
let oc =
open_out temp_file
in
get_input ()
|> output_string oc;
close_out oc;
Printf.sprintf "tar -zxf %s" temp_file
|> Sys.command |> ignore ;
let oc =
Filename.concat (Qpackage.ezfio_work ezfio_filename) "qp_run_address"
|> open_out
in
Printf.fprintf oc "tcp://%s:%d\n" localhost localport;
close_out oc;
Printf.printf " ...done\n%!"
end
end;
(* Main loop *)
let pollitem =
Zmq.Poll.mask_of [| (socket_in, Zmq.Poll.In) |]
in
let action () =
match Zmq.Socket.recv socket_in with
| "get_input" -> get_input ()
|> Zmq.Socket.send socket_in
| "get_ezfio_filename" -> get_ezfio_filename ()
|> 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
in
Printf.printf "
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)
(port+9) localhost (localport+9)
(Unix.gethostname ())
(port ) (localport )
(port+1) (localport+1)
(port+2) (localport+2)
(port+9) (localport+9);
Printf.printf "Ready\n%!";
while !run_status do
let polling =
Zmq.Poll.poll ~timeout:1000 pollitem
in
match polling.(0) with
| Some Zmq.Poll.In -> action ()
| None -> ()
| Some Zmq.Poll.In_out
| Some Zmq.Poll.Out -> ()
done;
let () =
match socket_out with
| Some socket_out -> Zmq.Socket.close socket_out
| None -> ()
in
Zmq.Socket.close socket_in
in
(Thread.create f) ()
in
(* Termination *)
Thread.join input_thread;
Thread.join fortran_thread;
Thread.join pub_thread;
Thread.join ocaml_thread;
Zmq.Context.terminate zmq_context;
Printf.printf "qp_tunnel exited properly.\n"

2
ocaml/qptypes_generator.ml
View File

@ -58,7 +58,7 @@ let input_data = "
* Det_number_max : int
assert (x > 0) ;
if (x > 50_00_000_000) then
if (x > 50_000_000_000) then
warning \"More than 50 billion determinants\";
* States_number : int

2
scripts/ezfio_interface/qp_edit_template
View File

@ -79,7 +79,7 @@ let get s =
| Ao_basis ->
f Ao_basis.(read, to_rst)
| Determinants_by_hand ->
f Determinants_by_hand.(read_maybe, to_rst)
f Determinants_by_hand.(read ~full:false, to_rst)
{section_to_rst}
end
with

7
src/ao_one_e_ints/README.rst
View File

@ -6,9 +6,8 @@ All the one-electron integrals in the |AO| basis are here.
The most important providers for usual quantum-chemistry calculation are:
* `ao_kinetic_integral` which are the kinetic operator integrals on the |AO| basis (see :file:`kin_ao_ints.irp.f`)
* `ao_nucl_elec_integral` which are the nuclear-elctron operator integrals on the |AO| basis (see :file:`pot_ao_ints.irp.f`)
* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis (see :file:`ao_mono_ints.irp.f`)
* `ao_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis
* `ao_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis
* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis
Note that you can find other interesting integrals related to the position operator in :file:`spread_dipole_ao.irp.f`.

41
src/bitmask/bitmask_cas_routines.irp.f
View File

@ -2,6 +2,29 @@ use bitmasks
integer function number_of_holes(key_in)
BEGIN_DOC
! Function that returns the number of holes in the inact space
!
! popcnt(
! xor(
! iand(
! reunion_of_core_inact_bitmask(1,1),
! xor(
! key_in(1,1),
! iand(
! key_in(1,1),
! cas_bitmask(1,1,1))
! )
! ),
! reunion_of_core_inact_bitmask(1,1)) )
!
! (key_in && cas_bitmask)
! +---------------------+
! electrons in cas xor key_in
! +---------------------------------+
! electrons outside of cas && reunion_of_core_inact_bitmask
! +------------------------------------------------------------------+
! electrons in the core/inact space xor reunion_of_core_inact_bitmask
! +---------------------------------------------------------------------------------+
! holes
END_DOC
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
@ -78,22 +101,6 @@ integer function number_of_holes(key_in)
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) )
else if(N_int == 8)then
number_of_holes = number_of_holes &
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(8,1), xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1)))), reunion_of_core_inact_bitmask(8,1)) )&
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(8,2), xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1)))), reunion_of_core_inact_bitmask(8,2)) )
else
do i = 1, N_int
number_of_holes = number_of_holes &
@ -108,7 +115,7 @@ integer function number_of_holes(key_in)
) &
) &
), reunion_of_core_inact_bitmask(i,1)) ) &
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) )
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1)))), reunion_of_core_inact_bitmask(i,2)) )
enddo
endif
end

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

@ -141,8 +141,8 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
b%cur=0
! ! Try to adjust n_tasks around nproc/2 seconds per job
! n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0)))
n_tasks = 1
n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0)))
! n_tasks = 1
end do
integer, external :: disconnect_from_taskserver

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

@ -85,7 +85,7 @@ subroutine run_selection_slave(thread,iproc,energy)
if(ctask > 0) then
call sort_selection_buffer(buf)
! call merge_selection_buffers(buf,buf2)
print *, task_id(1), pt2(1), buf%cur, ctask
!print *, task_id(1), pt2(1), buf%cur, ctask
call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask)
! buf%mini = buf2%mini
pt2(:) = 0d0

18
src/determinants/EZFIO.cfg
View File

@ -66,6 +66,12 @@ interface: ezfio
doc: Number of determinants in the current wave function
type: integer
[n_det_qp_edit]
interface: ezfio
doc: Number of determinants to print in qp_edit
type: integer
interface: ezfio
[psi_coef]
interface: ezfio
doc: Coefficients of the wave function
@ -78,6 +84,18 @@ doc: Determinants of the variational space
type: integer*8
size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det)
[psi_coef_qp_edit]
interface: ezfio
doc: Coefficients of the wave function
type: double precision
size: (determinants.n_det_qp_edit,determinants.n_states)
[psi_det_qp_edit]
interface: ezfio
doc: Determinants of the variational space
type: integer*8
size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det_qp_edit)
[expected_s2]
interface: ezfio
doc: Expected value of |S^2|

60
src/determinants/create_excitations.irp.f
View File

@ -12,6 +12,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
integer(bit_kind), intent(inout) :: key_in(N_int,2)
integer, intent(out) :: i_ok
integer :: k,j,i
integer(bit_kind) :: mask
use bitmasks
ASSERT (i_hole > 0 )
ASSERT (i_particle <= mo_num)
@ -19,31 +20,66 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
! hole
k = shiftr(i_hole-1,bit_kind_shift)+1
j = i_hole-shiftl(k-1,bit_kind_shift)-1
mask = ibset(0_bit_kind,j)
! check whether position j is occupied
if (ibits(key_in(k,ispin),j,1).eq.1) then
if (iand(key_in(k,ispin),mask) /= 0_bit_kind) then
key_in(k,ispin) = ibclr(key_in(k,ispin),j)
else
i_ok= -1
return
end if
! particle
k = shiftr(i_particle-1,bit_kind_shift)+1
j = i_particle-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin) = ibset(key_in(k,ispin),j)
mask = ibset(0_bit_kind,j)
if (iand(key_in(k,ispin),mask) == 0_bit_kind) then
key_in(k,ispin) = ibset(key_in(k,ispin),j)
else
i_ok= -1
return
end if
integer :: n_elec_tmp
n_elec_tmp = 0
do i = 1, N_int
n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
enddo
if(n_elec_tmp .ne. elec_num)then
!print*, n_elec_tmp,elec_num
!call debug_det(key_in,N_int)
i_ok = -1
endif
! integer :: n_elec_tmp
! n_elec_tmp = 0
! do i = 1, N_int
! n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
! enddo
! if(n_elec_tmp .ne. elec_num)then
! print*, n_elec_tmp,elec_num
! call debug_det(key_in,N_int)
! stop -1
! endif
end
subroutine build_singly_excited_wavefunction(i_hole,i_particle,ispin,det_out,coef_out)
implicit none
BEGIN_DOC
! Applies the single excitation operator : a^{dager}_(i_particle) a_(i_hole) of
! spin = ispin to the current wave function (psi_det, psi_coef)
END_DOC
integer, intent(in) :: i_hole,i_particle,ispin
integer(bit_kind), intent(out) :: det_out(N_int,2,N_det)
double precision, intent(out) :: coef_out(N_det,N_states)
integer :: k
integer :: i_ok
double precision :: phase
do k=1,N_det
coef_out(k,:) = psi_coef(k,:)
det_out(:,:,k) = psi_det(:,:,k)
call do_single_excitation(det_out(1,1,k),i_hole,i_particle,ispin,i_ok)
if (i_ok == 1) then
call get_phase(psi_det(1,1,k), det_out(1,1,k),phase,N_int)
coef_out(k,:) = phase * coef_out(k,:)
else
coef_out(k,:) = 0.d0
det_out(:,:,k) = psi_det(:,:,k)
endif
enddo
end
logical function is_spin_flip_possible(key_in,i_flip,ispin)
implicit none
BEGIN_DOC

1
src/determinants/density_matrix.irp.f
View File

@ -269,7 +269,6 @@ subroutine save_natural_mos
END_DOC
call set_natural_mos
call save_mos
end

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

@ -44,6 +44,16 @@ BEGIN_PROVIDER [ integer, N_det ]
ASSERT (N_det > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_qp_edit ]
implicit none
BEGIN_DOC
! Number of determinants to print in qp_edit
END_DOC
N_det_qp_edit = min(N_det,10000)
END_PROVIDER
BEGIN_PROVIDER [integer, max_degree_exc]
implicit none
integer :: i,degree
@ -476,7 +486,7 @@ subroutine save_wavefunction
endif
if (mpi_master) then
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
endif
endif
end
@ -504,12 +514,16 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
integer*8, allocatable :: psi_det_save(:,:,:)
double precision, allocatable :: psi_coef_save(:,:)
integer :: i,j,k
double precision :: accu_norm
integer :: i,j,k, ndet_qp_edit
if (mpi_master) then
ndet_qp_edit = min(ndet,N_det_qp_edit)
call ezfio_set_determinants_N_int(N_int)
call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_N_det(ndet)
call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
call ezfio_set_determinants_n_states(nstates)
call ezfio_set_determinants_mo_label(mo_label)
@ -522,10 +536,10 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
enddo
enddo
call ezfio_set_determinants_psi_det(psi_det_save)
call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
deallocate (psi_det_save)
allocate (psi_coef_save(ndet,nstates))
double precision :: accu_norm
do k=1,nstates
do i=1,ndet
psi_coef_save(i,k) = psicoef(i,k)
@ -535,6 +549,18 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
call ezfio_set_determinants_psi_coef(psi_coef_save)
deallocate (psi_coef_save)
allocate (psi_coef_save(ndet_qp_edit,nstates))
do k=1,nstates
do i=1,ndet_qp_edit
psi_coef_save(i,k) = psicoef(i,k)
enddo
call normalize(psi_coef_save(1,k),ndet_qp_edit)
enddo
call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save)
deallocate (psi_coef_save)
call write_int(6,ndet,'Saved determinants')
endif
end
@ -559,54 +585,80 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde
integer :: N_int2
equivalence (det_8, det_bk)
integer :: i,k
integer :: i,j,k, ndet_qp_edit
call ezfio_set_determinants_N_int(N_int)
call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_N_det(ndetsave)
call ezfio_set_determinants_n_states(nstates)
call ezfio_set_determinants_mo_label(mo_label)
if (mpi_master) then
ndet_qp_edit = min(ndetsave,N_det_qp_edit)
call ezfio_set_determinants_N_int(N_int)
call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_N_det(ndetsave)
call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
call ezfio_set_determinants_n_states(nstates)
call ezfio_set_determinants_mo_label(mo_label)
N_int2 = (N_int*bit_kind)/8
allocate (psi_det_save(N_int2,2,ndetsave))
do i=1,ndetsave
do k=1,N_int
det_bk(k) = psidet(k,1,index_det_save(i))
enddo
do k=1,N_int2
psi_det_save(k,1,i) = det_8(k)
enddo
do k=1,N_int
det_bk(k) = psidet(k,2,index_det_save(i))
enddo
do k=1,N_int2
psi_det_save(k,2,i) = det_8(k)
enddo
enddo
call ezfio_set_determinants_psi_det(psi_det_save)
deallocate (psi_det_save)
allocate (psi_coef_save(ndetsave,nstates))
double precision :: accu_norm(nstates)
accu_norm = 0.d0
do k=1,nstates
N_int2 = (N_int*bit_kind)/8
allocate (psi_det_save(N_int2,2,ndetsave))
do i=1,ndetsave
accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
psi_coef_save(i,k) = psicoef(index_det_save(i),k)
do k=1,N_int
det_bk(k) = psidet(k,1,index_det_save(i))
enddo
do k=1,N_int2
psi_det_save(k,1,i) = det_8(k)
enddo
do k=1,N_int
det_bk(k) = psidet(k,2,index_det_save(i))
enddo
do k=1,N_int2
psi_det_save(k,2,i) = det_8(k)
enddo
enddo
enddo
do k = 1, nstates
accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
enddo
do k=1,nstates
do i=1,ndetsave
psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
enddo
enddo
call ezfio_set_determinants_psi_det(psi_det_save)
call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
deallocate (psi_det_save)
call ezfio_set_determinants_psi_coef(psi_coef_save)
call write_int(6,ndet,'Saved determinants')
deallocate (psi_coef_save)
allocate (psi_coef_save(ndetsave,nstates))
double precision :: accu_norm(nstates)
accu_norm = 0.d0
do k=1,nstates
do i=1,ndetsave
accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
psi_coef_save(i,k) = psicoef(index_det_save(i),k)
enddo
enddo
do k = 1, nstates
accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
enddo
do k=1,nstates
do i=1,ndetsave
psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
enddo
enddo
call ezfio_set_determinants_psi_coef(psi_coef_save)
deallocate (psi_coef_save)
allocate (psi_coef_save(ndet_qp_edit,nstates))
accu_norm = 0.d0
do k=1,nstates
do i=1,ndet_qp_edit
accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
psi_coef_save(i,k) = psicoef(index_det_save(i),k)
enddo
enddo
do k = 1, nstates
accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
enddo
do k=1,nstates
do i=1,ndet_qp_edit
psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
enddo
enddo
call ezfio_set_determinants_psi_coef(psi_coef_save)
deallocate (psi_coef_save)
call write_int(6,ndet,'Saved determinants')
endif
end

12
src/determinants/occ_pattern.irp.f
View File

@ -53,7 +53,17 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Generate all possible determinants for a give occ_pattern
! Generate all possible determinants for a given occ_pattern
!
! Input :
! o : occupation pattern : (doubly occupied, singly occupied)
! sze : Number of produced determinants, computed by `occ_pattern_to_dets_size`
! n_alpha : Number of $\alpha$ electrons
! Nint : N_int
!
! Output:
! d : determinants
!
END_DOC
integer ,intent(in) :: Nint
integer ,intent(in) :: n_alpha ! Number of alpha electrons

15
src/dft_utils_in_r/kin_dens.irp.f Normal file
View File

@ -0,0 +1,15 @@
BEGIN_PROVIDER [double precision, kinetic_density_generalized, (n_points_final_grid)]
implicit none
integer :: i,j,m,i_point
kinetic_density_generalized = 0.d0
do i_point = 1, n_points_final_grid
do i = 1, mo_num
do j = 1, mo_num
do m = 1, 3
kinetic_density_generalized(i_point) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i_point) * mos_grad_in_r_array_tranp(m,i,i_point) * one_e_dm_mo_for_dft(j,i,1)
enddo
enddo
enddo
enddo
END_PROVIDER

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

@ -22,6 +22,7 @@
BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_num,n_points_final_grid,3)]
&BEGIN_PROVIDER[double precision, mos_grad_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! mos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith mo on the jth grid point
@ -35,8 +36,37 @@
do m=1,3
call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_num)
enddo
integer :: i,j
do i = 1, n_points_final_grid
do j = 1, mo_num
do m = 1, 3
mos_grad_in_r_array_tranp(m,j,i) = mos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)]
&BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)]
implicit none
integer :: i,m,j
alpha_dens_kin_in_r = 0.d0
beta_dens_kin_in_r = 0.d0
do i = 1, n_points_final_grid
do j = 1, elec_alpha_num
do m = 1, 3
alpha_dens_kin_in_r(i) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0
enddo
enddo
do j = 1, elec_beta_num
do m = 1, 3
beta_dens_kin_in_r(i) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_num,n_points_final_grid,3)]
implicit none
BEGIN_DOC

36
src/dft_utils_one_e/ec_lyp.irp.f
View File

@ -35,11 +35,8 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
include 'constants.include.F'
! Input variables
double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2
! Local variables
double precision :: a,b,c,d,c_f,omega,delta
double precision :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom
double precision :: thr,huge_num,rho_inv
@ -47,8 +44,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
double precision :: tmp1,tmp2,tmp3,tmp4
double precision :: big1,big2,big3
! Output variables
! Constants of the LYP correlation functional
@ -57,15 +52,25 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
c = 0.2533d0
d = 0.349d0
thr = 1d-10
ec_lyp_88 = 0.d0
thr = 1d-15
huge_num = 1.d0/thr
if(dabs(rho_a).lt.thr)then
return
endif
if(dabs(rho_b).lt.thr)then
return
endif
if(rho.lt.0.d0)then
print*,'pb !! rho.lt.0.d0'
stop
endif
rho_13 = rho**(1d0/3d0)
rho_113 = rho**(11d0/3d0)
rho_13 = rho**(1.d0/3.d0)
rho_113 = rho**(11.d0/3.d0)
if(dabs(rho_13) < thr) then
rho_inv_13 = huge_num
@ -76,13 +81,13 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
if (dabs(rho_113) < thr) then
rho_inv_113 = huge_num
else
rho_inv_113 = 1d0/rho_113
rho_inv_113 = 1.d0/rho_113
endif
if (dabs(rho) < thr) then
rho_inv = huge_num
else
rho_inv = 1d0/rho
rho_inv = 1.d0/rho
endif
! Useful quantities to predefine
@ -90,21 +95,21 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
denom = 1d0/(1d0 + d*rho_inv_13)
omega = rho_inv_113*exp(-c*rho_inv_13)*denom
delta = c*rho_inv_13 + d*rho_inv_13*denom
c_f = 0.3d0*(3d0*pi*pi)**(2d0/3d0)
c_f = 0.3d0*(3.d0*pi*pi)**(2.d0/3.d0)
rho_2 = rho *rho
rho_a_2 = rho_a*rho_a
rho_b_2 = rho_b*rho_b
cst_2_113 = 2d0**(11d0/3d0)
cst_8_3 = 8d0/3d0
cst_2_113 = 2.d0**(11.d0/3.d0)
cst_8_3 = 8.d0/3.d0
! first term in the equation (2) of Preuss CPL, 1989
big1 = 4d0*denom*rho_a*rho_b*rho_inv
big1 = 4.d0*denom*rho_a*rho_b*rho_inv
tmp1 = cst_2_113*c_f*(rho_a**cst_8_3 + rho_b**cst_8_3)
tmp2 = (47d0/18d0 - 7d0/18d0*delta)*grad_rho_2
tmp2 = (47.d0/18.d0 - 7.d0/18.d0*delta)*grad_rho_2
tmp3 = - (5d0/2d0 - 1.d0/18d0*delta)*(grad_rho_a_2 + grad_rho_b_2)
tmp4 = - (delta - 11d0)/9d0*(rho_a*rho_inv*grad_rho_a_2 + rho_b*rho_inv*grad_rho_b_2)
big2 = rho_a*rho_b*(tmp1 + tmp2 + tmp3 + tmp4)
@ -114,7 +119,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
tmp3 = grad_rho_a_2*(2d0/3d0*rho_2 - rho_b_2)
big3 = tmp1 + tmp2 + tmp3
ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3
end

100
src/dft_utils_one_e/ec_scan.irp.f Normal file
View File

@ -0,0 +1,100 @@
double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2)
include 'constants.include.F'
implicit none
double precision, intent(in) :: rho_a,rho_b,tau,grad_rho_2
double precision :: cst_13,cst_23,cst_43,cst_53,rho_inv,cst_18,cst_3pi2
double precision :: thr,nup,ndo,xi,s,spin_d,drho,drho2,rho,inv_1alph,e_c_lsda1,h0
double precision :: rs,t_w,t_unif,ds_xi,alpha,fc_alpha,step_f,cst_1alph,beta_inf
double precision :: c_1c,c_2c,d_c,e_c_ldsa1,h1,phi,t,beta_rs,gama,a,w_1,g_at2,phi_3,e_c_1
double precision :: b_1c,b_2c,b_3c,dx_xi,gc_xi,e_c_lsda0,w_0,g_inf,cx_xi,x_inf,f0,e_c_0
thr = 1.d-12
nup = max(rho_a,thr)
ndo = max(rho_b,thr)
rho = nup + ndo
ec_scan = 0.d0
if((rho).lt.thr)return
! constants ...
rho_inv = 1.d0/rho
cst_13 = 1.d0/3.d0
cst_23 = 2.d0 * cst_13
cst_43 = 4.d0 * cst_13
cst_53 = 5.d0 * cst_13
cst_18 = 1.d0/8.d0
cst_3pi2 = 3.d0 * pi*pi
drho2 = max(grad_rho_2,thr)
drho = dsqrt(drho2)
if((nup-ndo).gt.0.d0)then
spin_d = max(nup-ndo,thr)
else
spin_d = min(nup-ndo,-thr)
endif
c_1c = 0.64d0
c_2c = 1.5d0
d_c = 0.7d0
b_1c = 0.0285764d0
b_2c = 0.0889d0
b_3c = 0.125541d0
gama = 0.031091d0
! correlation energy lsda1
call ec_only_lda_sr(0.d0,nup,ndo,e_c_lsda1)
xi = spin_d/rho
rs = (cst_43 * pi * rho)**(-cst_13)
s = drho/( 2.d0 * cst_3pi2**(cst_13) * rho**cst_43 )
t_w = drho2 * cst_18 * rho_inv
ds_xi = 0.5d0 * ( (1.d0+xi)**cst_53 + (1.d0 - xi)**cst_53)
t_unif = 0.3d0 * (cst_3pi2)**cst_23 * rho**cst_53*ds_xi
t_unif = max(t_unif,thr)
alpha = (tau - t_w)/t_unif
cst_1alph= 1.d0 - alpha
if(cst_1alph.gt.0.d0)then
cst_1alph= max(cst_1alph,thr)
else
cst_1alph= min(cst_1alph,-thr)
endif
inv_1alph= 1.d0/cst_1alph
phi = 0.5d0 * ( (1.d0+xi)**cst_23 + (1.d0 - xi)**cst_23)
phi_3 = phi*phi*phi
t = (cst_3pi2/16.d0)**cst_13 * s / (phi * rs**0.5d0)
w_1 = dexp(-e_c_lsda1/(gama * phi_3)) - 1.d0
a = beta_rs(rs) /(gama * w_1)
g_at2 = 1.d0/(1.d0 + 4.d0 * a*t*t)**0.25d0
h1 = gama * phi_3 * dlog(1.d0 + w_1 * (1.d0 - g_at2))
! interpolation function
fc_alpha = dexp(-c_1c * alpha * inv_1alph) * step_f(cst_1alph) - d_c * dexp(c_2c * inv_1alph) * step_f(-cst_1alph)
! first part of the correlation energy
e_c_1 = e_c_lsda1 + h1
dx_xi = 0.5d0 * ( (1.d0+xi)**cst_43 + (1.d0 - xi)**cst_43)
gc_xi = (1.d0 - 2.3631d0 * (dx_xi - 1.d0) ) * (1.d0 - xi**12.d0)
e_c_lsda0= - b_1c / (1.d0 + b_2c * rs**0.5d0 + b_3c * rs)
w_0 = dexp(-e_c_lsda0/b_1c) - 1.d0
beta_inf = 0.066725d0 * 0.1d0 / 0.1778d0
cx_xi = -3.d0/(4.d0*pi) * (9.d0 * pi/4.d0)**cst_13 * dx_xi
x_inf = 0.128026d0
f0 = -0.9d0
g_inf = 1.d0/(1.d0 + 4.d0 * x_inf * s*s)**0.25d0
h0 = b_1c * dlog(1.d0 + w_0 * (1.d0 - g_inf))
e_c_0 = (e_c_lsda0 + h0) * gc_xi
ec_scan = e_c_1 + fc_alpha * (e_c_0 - e_c_1)
end
double precision function step_f(x)
implicit none
double precision, intent(in) :: x
if(x.lt.0.d0)then
step_f = 0.d0
else
step_f = 1.d0
endif
end
double precision function beta_rs(rs)
implicit none
double precision, intent(in) ::rs
beta_rs = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
end

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

@ -36,24 +36,24 @@ function run_stoch() {
[[ -n $TRAVIS ]] && skip
qp set_file f2.ezfio
qp set_frozen_core
run_stoch -199.30496 1.e-4
run_stoch -199.30486 1.e-4
}
@test "NH3" { # 10.6657s
qp set_file nh3.ezfio
qp set_mo_class --core="[1-4]" --act="[5-72]"
run -56.244753429144986 1.e-5
run -56.244753429144986 1.e-4
}
@test "DHNO" { # 11.4721s
qp set_file dhno.ezfio
qp set_mo_class --core="[1-7]" --act="[8-64]"
run -130.459020029816 1.e-5
run -130.459020029816 1.e-4
}
@test "HCO" { # 12.2868s
qp set_file hco.ezfio
run -113.297494345682 2.e-05
run -113.297494345682 1.e-4
}
@test "H2O2" { # 12.9214s
@ -65,82 +65,82 @@ function run_stoch() {
@test "HBO" { # 13.3144s
[[ -n $TRAVIS ]] && skip
qp set_file hbo.ezfio
run -100.214185815312 1.e-5
run -100.212829869715 1.e-4
}
@test "H2O" { # 11.3727s
[[ -n $TRAVIS ]] && skip
qp set_file h2o.ezfio
run -76.2359268957699 2.e-5
run -76.2359268957699 1.e-4
}
@test "ClO" { # 13.3755s
[[ -n $TRAVIS ]] && skip
qp set_file clo.ezfio
run -534.546005867797 5.e-5
run -534.545881614967 1.e-4
}
@test "SO" { # 13.4952s
[[ -n $TRAVIS ]] && skip
qp set_file so.ezfio
run -26.0124797722154 1.e-5
run -26.0158153138924 1.e-4
}
@test "H2S" { # 13.6745s
[[ -n $TRAVIS ]] && skip
qp set_file h2s.ezfio
run -398.859480581924 1.e-5
run -398.859168655255 1.e-4
}
@test "OH" { # 13.865s
[[ -n $TRAVIS ]] && skip
qp set_file oh.ezfio
run -75.6119887538831 1.e-05
run -75.6120779012574 1.e-4
}
@test "SiH2_3B1" { # 13.938ss
[[ -n $TRAVIS ]] && skip
qp set_file sih2_3b1.ezfio
run -290.017539006762 1.e-5
run -290.017539006762 1.e-4
}
@test "H3COH" { # 14.7299s
[[ -n $TRAVIS ]] && skip
qp set_file h3coh.ezfio
run -115.205054063687 1.e-5
run -115.205941463667 1.e-4
}
@test "SiH3" { # 15.99s
[[ -n $TRAVIS ]] && skip
qp set_file sih3.ezfio
run -5.57269434557089 2.e-05
run -5.57241217753818 1.e-4
}
@test "CH4" { # 16.1612s
[[ -n $TRAVIS ]] && skip
qp set_file ch4.ezfio
qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
run -40.2409059687324 2.e-5
run -40.2409678239136 1.e-4
}
@test "ClF" { # 16.8864s
[[ -n $TRAVIS ]] && skip
qp set_file clf.ezfio
run -559.170406471496 1.e-5
run -559.170272077166 1.e-4
}
@test "SO2" { # 17.5645s
[[ -n $TRAVIS ]] && skip
qp set_file so2.ezfio
qp set_mo_class --core="[1-8]" --act="[9-87]"
run -41.5746738713298 5.e-5
run -41.5746738713298 1.e-4
}
@test "C2H2" { # 17.6827s
[[ -n $TRAVIS ]] && skip
qp set_file c2h2.ezfio
qp set_mo_class --act="[1-30]" --del="[31-36]"
run -12.3670840202635 2.e-5
run -12.3656179738175 1.e-4
}
@test "N2" { # 18.0198s
@ -154,14 +154,14 @@ function run_stoch() {
[[ -n $TRAVIS ]] && skip
qp set_file n2h4.ezfio
qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
run -111.367234092521 2.e-5
run -111.367332681559 1.e-4
}
@test "CO2" { # 21.1748s
[[ -n $TRAVIS ]] && skip
qp set_file co2.ezfio
qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
run -187.969676381867 1.e-5
run -187.968599504402 1.e-4
}
@ -169,13 +169,13 @@ function run_stoch() {
[[ -n $TRAVIS ]] && skip
qp set_file cu_nh3_4_2plus.ezfio
qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
run -1862.98610987882 1.e-05
run -1862.98614665139 1.e-04
}
@test "HCN" { # 20.3273s
[[ -n $TRAVIS ]] && skip
qp set_file hcn.ezfio
qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
run -93.0799328685679 2.e-5
run -93.0728641601823 1.e-4
}

2
src/nuclei/EZFIO.cfg
View File

@ -21,7 +21,7 @@ type: double precision
size: (nuclei.nucl_num,3)
interface: ezfio
[disk_access_nuclear_repulsion]
[io_nuclear_repulsion]
doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ]
type: Disk_access
interface: ezfio,provider,ocaml

4
src/nuclei/nuclei.irp.f
View File

@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
END_DOC
PROVIDE mpi_master nucl_coord nucl_charge nucl_num
if (disk_access_nuclear_repulsion.EQ.'Read') then
if (io_nuclear_repulsion == 'Read') then
logical :: has
if (mpi_master) then
@ -194,7 +194,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
call write_time(6)
call write_double(6,nuclear_repulsion,'Nuclear repulsion energy')
if (disk_access_nuclear_repulsion.EQ.'Write') then
if (io_nuclear_repulsion == 'Write') then
if (mpi_master) then
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
endif

1
src/tools/print_ci_vectors.irp.f
View File

@ -24,6 +24,7 @@ subroutine routine
implicit none
integer :: i,k
integer :: degree
call print_energy_components
do i = 1, N_det
print *, 'Determinant ', i
call debug_det(psi_det(1,1,i),N_int)