LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-11 05:43:43 +02:00
Code Issues Releases Wiki Activity

Merge pull request #8 from LCPQ/master

Browse Source 
merge with LCPQ
This commit is contained in:
garniron 2016-02-17 12:27:34 +01:00
commit c41d5220c3
65 changed files with 4106 additions and 659 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.travis.yml
View File

@ -24,7 +24,7 @@ python:
script:
- ./configure --production ./config/gfortran.cfg
- source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD
- source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
- source ./quantum_package.rc ; ninja
- source ./quantum_package.rc ; cd ocaml ; make ; cd -
- source ./quantum_package.rc ; cd tests ; bats bats/qp.bats

33
README.md
View File

@ -11,6 +11,12 @@ Set of quantum chemistry programs and libraries.
For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions.
Demo
====
[![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo")
# Installation
@ -56,10 +62,14 @@ This file contains all the environment variables needed by the quantum package b
### Optional) Add some modules
Usage: qp_module.py list (--installed|--avalaible-local|--avalaible-remote)
qp_module.py install <name>...
qp_module.py create -n <name> [<children_module>...]
```
Usage:
qp_module.py create -n <name> [<children_modules>...]
qp_module.py download -n <name> [<path_folder>...]
qp_module.py install <name>...
qp_module.py list (--installed | --available-local)
qp_module.py uninstall <name>
```
For exemple you can type :
`qp_module.py install Full_CI`
@ -138,3 +148,20 @@ You have two or more ezfio configuration files for the same variable. Check file
- rm $QP_ROOT/install/EZFIO/config/*
- ninja
### Error: Seg Fault (139)
```
Segmentation fault (core dumped)
Program exited with code 139.
```
#### Why ?
It's caused when we call the DGEM routine of LAPACK.
##### Fix
Set `ulimit -s unlimited`, before runing `qp_run`. It seem to fix the problem.

4
configure vendored
View File

@ -533,13 +533,13 @@ def recommendation():
print " source {0}".format(path)
print ""
print "Then, install the modules you want to install using :"
print " qp_install_module.py "
print " qp_module.py"
print ""
print "Finally :"
print " ninja"
print " make -C ocaml"
print ""
print "You can install more plugin with the qp_install command"
print "You can install more plugin with the qp_module.py install command"
print "PS : For more info on compiling the code, read the README.md"

168
data/basis/chipman-dzp Normal file
View File

@ -0,0 +1,168 @@
HYDROGEN
S 3
1 127.9500000 0.0107360
2 19.2406000 0.1195440
3 2.8992000 0.9264160
S 1
1 0.6534000 1.0000000
S 1
1 0.1776000 1.0000000
S 1
1 0.0483000 1.0000000
P 1
1 1.0000000 1.0000000
BORON
S 5
1 2788.4100000 0.0021220
2 419.0390000 0.0161710
3 96.4683000 0.0783560
4 28.0694000 0.2632500
5 9.3760000 0.5967290
S 1
1 1.3057000 1.0000000
S 1
1 3.4062000 1.0000000
S 1
1 0.3245000 1.0000000
S 1
1 0.1022000 1.0000000
S 1
1 0.0330000 1.0000000
P 4
1 11.3413000 0.0179870
2 2.4360000 0.1103390
3 0.6836000 0.3831110
4 0.2134000 0.6478600
P 1
1 0.0701000 1.0000000
P 1
1 0.0226000 1.0000000
D 1
1 0.1600000 1.0000000
D 1
1 0.6400000 1.0000000
CARBON
S 5
1 5909.4400000 0.0020040
2 887.4510000 0.0153100
3 204.7490000 0.0742930
4 59.8376000 0.2533640
5 19.9981000 0.6005760
S 1
1 2.6860000 1.0000000
S 1
1 7.1927000 1.0000000
S 1
1 0.7000000 1.0000000
S 1
1 0.2133000 1.0000000
S 1
1 0.0667000 1.0000000
P 4
1 26.7860000 0.0182570
2 5.9564000 0.1164070
3 1.7074000 0.3901110
4 0.5314000 0.6372210
P 1
1 0.1654000 1.0000000
P 1
1 0.0517000 1.0000000
D 1
1 0.3700000 1.0000000
D 1
1 1.4800000 1.0000000
NITROGEN
S 5
1 5909.4400000 0.0020040
2 887.4510000 0.0153100
3 204.7490000 0.0742930
4 59.8376000 0.2533640
5 19.9981000 0.6005760
S 1
1 2.6860000 1.0000000
S 1
1 7.1927000 1.0000000
S 1
1 0.7000000 1.0000000
S 1
1 0.2133000 1.0000000
S 1
1 0.0667000 1.0000000
P 4
1 26.7860000 0.0182570
2 5.9564000 0.1164070
3 1.7074000 0.3901110
4 0.5314000 0.6372210
P 1
1 0.1654000 1.0000000
P 1
1 0.0517000 1.0000000
D 1
1 0.3700000 1.0000000
D 1
1 1.4800000 1.0000000
OXYGEN
S 5
1 7816.5400000 0.0020310
2 1175.8200000 0.0154360
3 273.1880000 0.0737710
4 81.1696000 0.2476060
5 27.1836000 0.6118320
S 1
1 3.4136000 1.0000000
S 1
1 9.5322000 1.0000000
S 1
1 0.9398000 1.0000000
S 1
1 0.2846000 1.0000000
S 1
1 0.0862000 1.0000000
P 4
1 35.1832000 0.0195800
2 7.9040000 0.1241890
3 2.3051000 0.3947270
4 0.7171000 0.6273750
P 1
1 0.2137000 1.0000000
P 1
1 0.0648000 1.0000000
D 1
1 0.5500000 1.0000000
D 1
1 2.2000000 1.0000000
FLUORINE
S 5
1 9994.7900000 0.0020170
2 1506.0300000 0.0152950
3 350.2690000 0.0731100
4 104.0530000 0.2464200
5 34.8432000 0.6125930
S 1
1 4.3688000 1.0000000
S 1
1 12.2164000 1.0000000
S 1
1 1.2078000 1.0000000
S 1
1 0.3634000 1.0000000
S 1
1 0.1101000 1.0000000
P 4
1 44.3555000 0.0208680
2 10.0820000 0.1300920
3 2.9959000 0.3962190
4 0.9383000 0.6203680
P 1
1 0.2733000 1.0000000
P 1
1 0.0828000 1.0000000
D 1
1 0.6650000 1.0000000
D 1
1 2.6600000 1.0000000

20
ocaml/Qputils.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core.Std
(*
let rec transpose = function
@ -24,5 +24,21 @@ let input_to_sexp s =
print_endline ("("^result^")");
"("^result^")"
|> Sexp.of_string
;;
let rmdir dirname =
let rec remove_one dir =
Sys.chdir dir;
Sys.readdir "."
|> Array.iter ~f:(fun x ->
match (Sys.is_directory x, Sys.is_file x) with
| (`Yes, _) -> remove_one x
| (_, `Yes) -> Sys.remove x
| _ -> failwith ("Unable to remove file "^x^".")
);
Sys.chdir "..";
Unix.rmdir dir
in
remove_one dirname

555
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -1,6 +1,6 @@
open Qputils;;
open Qptypes;;
open Core.Std;;
open Qputils
open Qptypes
open Core.Std
let spec =
let open Command.Spec in
@ -85,7 +85,7 @@ let list_basis () =
)
in
List.sort basis_list ~cmp:String.ascending
|> String.concat ~sep:"\t"
|> String.concat ~sep:"\n"
(** Run the program *)
@ -316,289 +316,316 @@ let run ?o b c d m p cart xyz_file =
if Sys.file_exists_exn ezfio_file then
failwith (ezfio_file^" already exists");
(* Create EZFIO *)
Ezfio.set_file ezfio_file;
let write_file () =
(* Create EZFIO *)
Ezfio.set_file ezfio_file;
(* Write Pseudo *)
let pseudo =
List.map nuclei ~f:(fun x ->
match pseudo_channel x.Atom.element with
| Some channel -> Pseudo.read_element channel x.Atom.element
| None -> Pseudo.empty x.Atom.element
)
in
(* Write Pseudo *)
let pseudo =
List.map nuclei ~f:(fun x ->
match pseudo_channel x.Atom.element with
| Some channel -> Pseudo.read_element channel x.Atom.element
| None -> Pseudo.empty x.Atom.element
)
in
let molecule =
let n_elec_to_remove =
List.fold pseudo ~init:0 ~f:(fun accu x ->
accu + (Positive_int.to_int x.Pseudo.n_elec))
in
{ Molecule.elec_alpha =
(Elec_alpha_number.to_int molecule.Molecule.elec_alpha)
- n_elec_to_remove/2
|> Elec_alpha_number.of_int;
Molecule.elec_beta =
(Elec_beta_number.to_int molecule.Molecule.elec_beta)
- (n_elec_to_remove - n_elec_to_remove/2)
|> Elec_beta_number.of_int;
Molecule.nuclei =
let charges =
List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec
|> Float.of_int)
|> Array.of_list
let molecule =
let n_elec_to_remove =
List.fold pseudo ~init:0 ~f:(fun accu x ->
accu + (Positive_int.to_int x.Pseudo.n_elec))
in
List.mapi molecule.Molecule.nuclei ~f:(fun i x ->
{ x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i)
|> Charge.of_float }
)
}
in
let nuclei =
molecule.Molecule.nuclei @ dummy
in
(* Write Electrons *)
Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int
molecule.Molecule.elec_alpha ) ;
Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int
molecule.Molecule.elec_beta ) ;
(* Write Nuclei *)
let labels =
List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei
and charges =
List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei
and coords =
(List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @
(List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @
(List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in
let nucl_num = (List.length labels) in
Ezfio.set_nuclei_nucl_num nucl_num ;
Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| nucl_num |] ~data:labels);
Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| nucl_num |] ~data:charges);
Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords);
(* Write pseudopotential *)
let () =
match p with
| None -> Ezfio.set_pseudo_do_pseudo false
| _ -> Ezfio.set_pseudo_do_pseudo true
in
let klocmax =
List.fold pseudo ~init:0 ~f:(fun accu x ->
let x =
List.length x.Pseudo.local
{ Molecule.elec_alpha =
(Elec_alpha_number.to_int molecule.Molecule.elec_alpha)
- n_elec_to_remove/2
|> Elec_alpha_number.of_int;
Molecule.elec_beta =
(Elec_beta_number.to_int molecule.Molecule.elec_beta)
- (n_elec_to_remove - n_elec_to_remove/2)
|> Elec_beta_number.of_int;
Molecule.nuclei =
let charges =
List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec
|> Float.of_int)
|> Array.of_list
in
List.mapi molecule.Molecule.nuclei ~f:(fun i x ->
{ x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i)
|> Charge.of_float }
)
}
in
if (x > accu) then x
else accu
)
and kmax =
List.fold pseudo ~init:0 ~f:(fun accu x ->
let x =
List.length x.Pseudo.non_local
let nuclei =
molecule.Molecule.nuclei @ dummy
in
if (x > accu) then x
else accu
)
and lmax =
List.fold pseudo ~init:0 ~f:(fun accu x ->
let x =
List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) ->
(* Write Electrons *)
Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int
molecule.Molecule.elec_alpha ) ;
Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int
molecule.Molecule.elec_beta ) ;
(* Write Nuclei *)
let labels =
List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei
and charges =
List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei
and coords =
(List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @
(List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @
(List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in
let nucl_num = (List.length labels) in
Ezfio.set_nuclei_nucl_num nucl_num ;
Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| nucl_num |] ~data:labels);
Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| nucl_num |] ~data:charges);
Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords);
(* Write pseudopotential *)
let () =
match p with
| None -> Ezfio.set_pseudo_do_pseudo false
| _ -> Ezfio.set_pseudo_do_pseudo true
in
let klocmax =
List.fold pseudo ~init:0 ~f:(fun accu x ->
let x =
Positive_int.to_int x.Pseudo.Primitive_non_local.proj
List.length x.Pseudo.local
in
if (x > accu) then x
else accu
)
and lmax =
List.fold pseudo ~init:0 ~f:(fun accu x ->
let x =
List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) ->
let x =
Positive_int.to_int x.Pseudo.Primitive_non_local.proj
in
if (x > accu) then x
else accu
)
in
if (x > accu) then x
else accu
)
in
if (x > accu) then x
else accu
)
in
let () =
Ezfio.set_pseudo_pseudo_klocmax klocmax;
Ezfio.set_pseudo_pseudo_kmax kmax;
Ezfio.set_pseudo_pseudo_lmax lmax;
let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k =
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0
in
List.iteri pseudo ~f:(fun j x ->
List.iteri x.Pseudo.local ~f:(fun i (y,c) ->
tmp_array_v_k.(i).(j) <- AO_coef.to_float c;
let y, z =
AO_expo.to_float y.Pseudo.Primitive_local.expo,
R_power.to_int y.Pseudo.Primitive_local.r_power
let kmax =
Array.init (lmax+1) ~f:(fun i->
List.map pseudo ~f:(fun x ->
List.filter x.Pseudo.non_local ~f:(fun (y,_) ->
(Positive_int.to_int y.Pseudo.Primitive_non_local.proj) = i)
|> List.length )
|> List.fold ~init:0 ~f:(fun accu x ->
if accu > x then accu else x)
)
|> Array.fold ~init:0 ~f:(fun accu i ->
if i > accu then i else accu)
in
let () =
Ezfio.set_pseudo_pseudo_klocmax klocmax;
Ezfio.set_pseudo_pseudo_kmax kmax;
Ezfio.set_pseudo_pseudo_lmax lmax;
let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k =
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0
in
tmp_array_dz_k.(i).(j) <- y;
tmp_array_n_k.(i).(j) <- z;
)
);
let concat_2d tmp_array =
let data =
Array.map tmp_array ~f:Array.to_list
|> Array.to_list
|> List.concat
in
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data
in
concat_2d tmp_array_v_k
|> Ezfio.set_pseudo_pseudo_v_k ;
concat_2d tmp_array_dz_k
|> Ezfio.set_pseudo_pseudo_dz_k;
concat_2d tmp_array_n_k
|> Ezfio.set_pseudo_pseudo_n_k;
List.iteri pseudo ~f:(fun j x ->
List.iteri x.Pseudo.local ~f:(fun i (y,c) ->
tmp_array_v_k.(i).(j) <- AO_coef.to_float c;
let y, z =
AO_expo.to_float y.Pseudo.Primitive_local.expo,
R_power.to_int y.Pseudo.Primitive_local.r_power
in
tmp_array_dz_k.(i).(j) <- y;
tmp_array_n_k.(i).(j) <- z;
)
);
let concat_2d tmp_array =
let data =
Array.map tmp_array ~f:Array.to_list
|> Array.to_list
|> List.concat
in
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data
in
concat_2d tmp_array_v_k
|> Ezfio.set_pseudo_pseudo_v_k ;
concat_2d tmp_array_dz_k
|> Ezfio.set_pseudo_pseudo_dz_k;
concat_2d tmp_array_n_k
|> Ezfio.set_pseudo_pseudo_n_k;
let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl =
Array.create ~len:(lmax+1)
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ),
Array.create ~len:(lmax+1)
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ),
Array.create ~len:(lmax+1)
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 )
in
List.iteri pseudo ~f:(fun j x ->
List.iteri x.Pseudo.non_local ~f:(fun i (y,c) ->
let k, y, z =
Positive_int.to_int y.Pseudo.Primitive_non_local.proj,
AO_expo.to_float y.Pseudo.Primitive_non_local.expo,
R_power.to_int y.Pseudo.Primitive_non_local.r_power
let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl =
Array.init (lmax+1) ~f:(fun _ ->
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )),
Array.init (lmax+1) ~f:(fun _ ->
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )),
Array.init (lmax+1) ~f:(fun _ ->
(Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 ))
in
tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c;
tmp_array_dz_kl.(k).(i).(j) <- y;
tmp_array_n_kl.(k).(i).(j) <- z;
)
);
let concat_3d tmp_array =
let data =
Array.map tmp_array ~f:(fun x ->
Array.map x ~f:Array.to_list
|> Array.to_list
|> List.concat)
|> Array.to_list
List.iteri pseudo ~f:(fun j x ->
let last_idx =
Array.create ~len:(lmax+1) 0
in
List.iter x.Pseudo.non_local ~f:(fun (y,c) ->
let k, y, z =
Positive_int.to_int y.Pseudo.Primitive_non_local.proj,
AO_expo.to_float y.Pseudo.Primitive_non_local.expo,
R_power.to_int y.Pseudo.Primitive_non_local.r_power
in
let i =
last_idx.(k)
in
tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c;
tmp_array_dz_kl.(k).(i).(j) <- y;
tmp_array_n_kl.(k).(i).(j) <- z;
last_idx.(k) <- i+1;
)
);
let concat_3d tmp_array =
let data =
Array.map tmp_array ~f:(fun x ->
Array.map x ~f:Array.to_list
|> Array.to_list
|> List.concat)
|> Array.to_list
|> List.concat
in
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data
in
concat_3d tmp_array_v_kl
|> Ezfio.set_pseudo_pseudo_v_kl ;
concat_3d tmp_array_dz_kl
|> Ezfio.set_pseudo_pseudo_dz_kl ;
concat_3d tmp_array_n_kl
|> Ezfio.set_pseudo_pseudo_n_kl ;
in
(* Write Basis set *)
let basis =
let nmax =
Nucl_number.get_max ()
in
let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
| [] -> accu
| e::tail ->
let new_accu =
(e,(Nucl_number.of_int ~max:nmax n))::accu
in
do_work new_accu (n+1) tail
in
let result = do_work [] 1 nuclei
|> List.rev
|> List.map ~f:(fun (x,i) ->
try
let e =
match x.Atom.element with
| Element.X -> Element.H
| e -> e
in
Basis.read_element (basis_channel x.Atom.element) i e
with
| End_of_file -> failwith
("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
)
|> List.concat
in
(* close all in_channels *)
result
in
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data
in
concat_3d tmp_array_v_kl
|> Ezfio.set_pseudo_pseudo_v_kl ;
concat_3d tmp_array_dz_kl
|> Ezfio.set_pseudo_pseudo_dz_kl ;
concat_3d tmp_array_n_kl
|> Ezfio.set_pseudo_pseudo_n_kl ;
in
(* Write Basis set *)
let basis =
let nmax =
Nucl_number.get_max ()
in
let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
| [] -> accu
| e::tail ->
let new_accu =
(e,(Nucl_number.of_int ~max:nmax n))::accu
let long_basis = Long_basis.of_basis basis in
let ao_num = List.length long_basis in
Ezfio.set_ao_basis_ao_num ao_num;
Ezfio.set_ao_basis_ao_basis b;
let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc)
and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n)
and ao_power=
let l = List.map long_basis ~f:(fun (x,_,_) -> x) in
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) )
in
do_work new_accu (n+1) tail
in
let result = do_work [] 1 nuclei
|> List.rev
|> List.map ~f:(fun (x,i) ->
try
let e =
match x.Atom.element with
| Element.X -> Element.H
| e -> e
in
Basis.read_element (basis_channel x.Atom.element) i e
with
| End_of_file -> failwith
("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
)
|> List.concat
in
(* close all in_channels *)
result
in
let long_basis = Long_basis.of_basis basis in
let ao_num = List.length long_basis in
Ezfio.set_ao_basis_ao_num ao_num;
Ezfio.set_ao_basis_ao_basis b;
let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc)
and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n)
and ao_power=
let l = List.map long_basis ~f:(fun (x,_,_) -> x) in
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@
(List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) )
in
let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x ->
if x > s then x
else s) ao_prim_num
in
let gtos =
List.map long_basis ~f:(fun (_,x,_) -> x)
in
let create_expo_coef ec =
let coefs =
begin match ec with
| `Coefs -> List.map gtos ~f:(fun x->
List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) )
| `Expos -> List.map gtos ~f:(fun x->
List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float
prim.Primitive.expo) )
end
let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x ->
if x > s then x
else s) ao_prim_num
in
let rec get_n n accu = function
| [] -> List.rev accu
| h::tail ->
let y =
begin match List.nth h n with
| Some x -> x
| None -> 0.
let gtos =
List.map long_basis ~f:(fun (_,x,_) -> x)
in
let create_expo_coef ec =
let coefs =
begin match ec with
| `Coefs -> List.map gtos ~f:(fun x->
List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) )
| `Expos -> List.map gtos ~f:(fun x->
List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float
prim.Primitive.expo) )
end
in
get_n n (y::accu) tail
in
let rec get_n n accu = function
| [] -> List.rev accu
| h::tail ->
let y =
begin match List.nth h n with
| Some x -> x
| None -> 0.
end
in
get_n n (y::accu) tail
in
let rec build accu = function
| n when n=ao_prim_num_max -> accu
| n -> build ( accu @ (get_n n [] coefs) ) (n+1)
in
build [] 0
in
let rec build accu = function
| n when n=ao_prim_num_max -> accu
| n -> build ( accu @ (get_n n [] coefs) ) (n+1)
in
build [] 0
in
let ao_coef = create_expo_coef `Coefs
and ao_expo = create_expo_coef `Expos
let ao_coef = create_expo_coef `Coefs
and ao_expo = create_expo_coef `Expos
in
let () =
Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ;
Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ;
Ezfio.set_ao_basis_ao_cartesian(cart);
in
match Input.Ao_basis.read () with
| None -> failwith "Error in basis"
| Some x -> Input.Ao_basis.write x
in
let () =
Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ;
Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ;
Ezfio.set_ao_basis_ao_cartesian(cart);
in
match Input.Ao_basis.read () with
| None -> failwith "Error in basis"
| Some x -> Input.Ao_basis.write x
try write_file () with
| ex ->
begin
begin
match Sys.is_directory ezfio_file with
| `Yes -> rmdir ezfio_file
| _ -> ()
end;
raise ex;
end
in ()

36
ocaml/qp_edit.ml
View File

@ -17,12 +17,12 @@ type keyword =
| Electrons
| Mo_basis
| Nuclei
| Determinants
| Integrals_bielec
| Hartree_fock
| Pseudo
| Integrals_bielec
| Perturbation
| Properties
| Hartree_fock
| Determinants
;;
@ -32,12 +32,12 @@ let keyword_to_string = function
| Electrons -> "Electrons"
| Mo_basis -> "MO basis"
| Nuclei -> "Molecule"
| Determinants -> "Determinants"
| Integrals_bielec -> "Integrals_bielec"
| Hartree_fock -> "Hartree_fock"
| Pseudo -> "Pseudo"
| Integrals_bielec -> "Integrals_bielec"
| Perturbation -> "Perturbation"
| Properties -> "Properties"
| Hartree_fock -> "Hartree_fock"
| Determinants -> "Determinants"
;;
@ -86,18 +86,18 @@ let get s =
f Ao_basis.(read, to_rst)
| Determinants_by_hand ->
f Determinants_by_hand.(read, to_rst)
| Determinants ->
f Determinants.(read, to_rst)
| Integrals_bielec ->
f Integrals_bielec.(read, to_rst)
| Hartree_fock ->
f Hartree_fock.(read, to_rst)
| Pseudo ->
f Pseudo.(read, to_rst)
| Integrals_bielec ->
f Integrals_bielec.(read, to_rst)
| Perturbation ->
f Perturbation.(read, to_rst)
| Properties ->
f Properties.(read, to_rst)
| Hartree_fock ->
f Hartree_fock.(read, to_rst)
| Determinants ->
f Determinants.(read, to_rst)
end
with
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
@ -135,12 +135,12 @@ let set str s =
in
let open Input in
match s with
| Determinants -> write Determinants.(of_rst, write) s
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
| Pseudo -> write Pseudo.(of_rst, write) s
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
| Perturbation -> write Perturbation.(of_rst, write) s
| Properties -> write Properties.(of_rst, write) s
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
| Determinants -> write Determinants.(of_rst, write) s
| Electrons -> write Electrons.(of_rst, write) s
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
| Nuclei -> write Nuclei.(of_rst, write) s
@ -188,12 +188,12 @@ let run check_only ezfio_filename =
Nuclei ;
Ao_basis;
Electrons ;
Determinants ;
Integrals_bielec ;
Hartree_fock ;
Pseudo ;
Integrals_bielec ;
Perturbation ;
Properties ;
Hartree_fock ;
Determinants ;
Mo_basis;
Determinants_by_hand ;
]

5
plugins/All_singles/EZFIO.cfg Normal file
View File

@ -0,0 +1,5 @@
[energy]
type: double precision
doc: Calculated Selected all_singles or all_1h_1p energy
interface: ezfio

18
plugins/All_singles/H_apply.irp.f Normal file
View File

@ -0,0 +1,18 @@
use bitmasks
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("just_1h_1p")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.filter_only_1h1p()
print s
s = H_apply("just_mono")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_double_excitations()
print s
END_SHELL

1
plugins/All_singles/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Generators_restart Perturbation Properties Selectors_no_sorted Utils

12
plugins/All_singles/README.rst Normal file
View File

@ -0,0 +1,12 @@
===========
All_singles
===========
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

76
plugins/All_singles/all_1h_1p.irp.f Normal file
View File

@ -0,0 +1,76 @@
program restart_more_singles
BEGIN_DOC
! Generates and select single and double excitations of type 1h-1p
! on the top of a given restart wave function of type CAS
END_DOC
read_wf = .true.
touch read_wf
print*,'ref_bitmask_energy = ',ref_bitmask_energy
call routine
end
subroutine routine
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:)
integer :: N_st, degree
integer :: n_det_before
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
i = 0
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
pt2=-1.d0
E_before = ref_bitmask_energy
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_1h_1p(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
E_before = CI_energy
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
call save_wavefunction
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))
enddo
endif
call ezfio_set_all_singles_energy(CI_energy)
call save_wavefunction
deallocate(pt2,norm_pert)
end

76
plugins/All_singles/all_singles.irp.f Normal file
View File

@ -0,0 +1,76 @@
program restart_more_singles
BEGIN_DOC
! Generates and select single excitations
! on the top of a given restart wave function
END_DOC
read_wf = .true.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision,allocatable :: E_before(:)
integer :: n_det_before
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
i = 0
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
pt2=-1.d0
E_before = ref_bitmask_energy
pt2_max = 1.d-10
n_det_max = 200000
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy
print*,'pt2 = ',pt2
print*,'E+PT2 = ',E_before + pt2
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
E_before = CI_energy
call save_wavefunction
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))
enddo
endif
call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end

16
plugins/Full_CI/full_ci.irp.f
View File

@ -62,11 +62,27 @@ program full_ci
endif
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k = 1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E(before)+PT2 = ', E_CI_before+pt2
enddo
print *, '-----'
E_CI_before = CI_energy
if(N_states.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_states
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_states
print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
enddo
endif
E_CI_before = CI_energy
call ezfio_set_full_ci_energy(CI_energy)
if (abort_all) then
exit

2
plugins/Hartree_Fock/.gitignore vendored
View File

@ -5,7 +5,6 @@ AO_Basis
Bitmask
Electrons
Ezfio_files
Huckel_guess
IRPF90_man
IRPF90_temp
Integrals_Bielec
@ -16,7 +15,6 @@ Makefile
Makefile.depend
Nuclei
Pseudo
SCF
Utils
ZMQ
ezfio_interface.irp.f

7
plugins/Hartree_Fock/EZFIO.cfg
View File

@ -26,3 +26,10 @@ default: Huckel
type: double precision
doc: Calculated HF energy
interface: ezfio
[no_oa_or_av_opt]
type: logical
doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the SCF procedure
interface: ezfio,provider,ocaml
default: False

2
plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Integrals_Bielec MOGuess
Integrals_Bielec MOGuess Bitmask

74
plugins/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -11,63 +11,35 @@
double precision, allocatable :: work(:), F(:,:), S(:,:)
! if (mo_tot_num == ao_num) then
! ! Solve H.C = E.S.C in AO basis set
!
! allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
! do j=1,ao_num
! do i=1,ao_num
! S(i,j) = ao_overlap(i,j)
! F(i,j) = Fock_matrix_ao(i,j)
! enddo
! enddo
!
! n = ao_num
! lwork = 1+6*n + 2*n*n
! liwork = 3 + 5*n
!
! allocate(work(lwork), iwork(liwork) )
!
! lwork = -1
! liwork = -1
!
! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
!
! if (info /= 0) then
! print *, irp_here//' failed : ', info
! stop 1
! endif
! lwork = int(work(1))
! liwork = iwork(1)
! deallocate(work,iwork)
! allocate(work(lwork), iwork(liwork) )
!
! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
!
! if (info /= 0) then
! print *, irp_here//' failed : ', info
! stop 1
! endif
! do j=1,mo_tot_num
! do i=1,ao_num
! eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
! enddo
! enddo
!
! deallocate(work, iwork, F, S)
!
! else
!
! Solve H.C = E.C in MO basis set
allocate( F(mo_tot_num_align,mo_tot_num) )
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
enddo
enddo
if(no_oa_or_av_opt)then
integer :: iorb,jorb
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
enddo
endif
! Insert level shift here

22
plugins/MRCC_Utils/H_apply.irp.f
View File

@ -24,5 +24,27 @@ s.data["size_max"] = "3072"
print s
s = H_apply("mrcepa")
s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
s.data["declarations"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
"""
s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
s.data["decls_main"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
"""
s.data["finalization"] = ""
s.data["copy_buffer"] = ""
s.data["generate_psi_guess"] = ""
s.data["size_max"] = "3072"
# print s
END_SHELL

4
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -26,7 +26,11 @@
! TODO --- Test perturbatif ------
do k=1,N_states
lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
! TODO : i_h_psi peut sortir de la boucle?
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
if (ihpsi_current(k) == 0.d0) then
ihpsi_current(k) = 1.d-32
endif
tmp = psi_non_ref_coef(i,k)/ihpsi_current(k)
i_pert = 0
! Perturbation only if 1st order < 0.5 x second order

260
plugins/MRCC_Utils/mrcepa_dress.irp.f Normal file
View File

@ -0,0 +1,260 @@
use omp_lib
use bitmasks
subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc
integer, intent(in) :: Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l
integer :: degree_alpha(psi_det_size)
integer :: idx_alpha(0:psi_det_size)
logical :: good, fullMatch
integer(bit_kind) :: tq(Nint,2,n_selected)
integer :: N_tq, c_ref ,degree
double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
double precision, allocatable :: dIa_hia(:,:)
double precision :: haj, phase, phase2
double precision :: f(N_states), ci_inv(N_states)
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(Nint,2)
integer(bit_kind) :: tmp_det_0(Nint,2)
integer :: iint, ipos
integer :: i_state, i_sd, k_sd, l_sd, i_I, i_alpha
integer(bit_kind),allocatable :: miniList(:,:,:)
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer,allocatable :: idx_miniList(:)
integer :: N_miniList, ni, leng
integer(bit_kind) :: isum
double precision :: hia
integer, allocatable :: index_sorted(:)
leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det))
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
if(fullMatch) then
return
end if
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
allocate (dIa_hia(N_states,Ndet_non_ref))
! |I>
! |alpha>
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
end if
do i_alpha=1,N_tq
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
integer, external :: get_index_in_psi_det_sorted_bit
index_sorted = huge(-1)
do j=1,idx_alpha(0)
idx_alpha(j) = idx_miniList(idx_alpha(j))
index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j)
end do
! |I>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
if (degree > 4) then
cycle
endif
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
!TODO: MR
do i_sd=1,idx_alpha(0)
call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint)
if (degree > 2) then
cycle
endif
call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
tmp_det_0 = 0_bit_kind
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
endif
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia)
! <I| <> |alpha>
do k_sd=1,idx_alpha(0)
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
if (degree > 2) then
cycle
endif
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
tmp_det = 0_bit_kind
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
isum = 0_bit_kind
do iint = 1,N_int
isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) &
+ iand(tmp_det(iint,2), tmp_det_0(iint,2))
enddo
if (isum /= 0_bit_kind) then
cycle
endif
! <I| /k\ |alpha>
! <I|H|k>
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
do i_state=1,N_states
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do k=1,N_int
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
! <I| \l/ |alpha>
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) )
! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl)
! do i_state=1,N_states
! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2
! enddo
do l_sd=1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
if (degree == 0) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
do i_state=1,N_states
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
enddo
exit
endif
enddo
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
enddo
do i_state=1,N_states
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
enddo
k_sd = idx_alpha(i_sd)
do i_state=1,N_states
dIa_hia(i_state,k_sd) = dIa(i_state) * hia
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do i_state=1,N_states
delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd)
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state)
else
delta_ii_(i_I,i_state) = 0.d0
endif
enddo
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
enddo
deallocate (dIa_hia,index_sorted)
deallocate(miniList, idx_miniList)
end

97
plugins/MRCC_Utils/mrcepa_general.irp.f Normal file
View File

@ -0,0 +1,97 @@
subroutine run_mrcepa
implicit none
call set_generators_bitmasks_as_holes_and_particles
call mrcepa_iterations
end
subroutine mrcepa_iterations
implicit none
integer :: i,j
double precision :: E_new, E_old, delta_e
integer :: iteration,i_oscillations
double precision :: E_past(4)
E_new = 0.d0
delta_E = 1.d0
iteration = 0
j = 1
i_oscillations = 0
do while (delta_E > 1.d-7)
iteration += 1
print *, '==========================='
print *, 'MRCEPA Iteration', iteration
print *, '==========================='
print *, ''
E_old = sum(ci_energy_dressed)
call write_double(6,ci_energy_dressed(1),"MRCEPA energy")
call diagonalize_ci_dressed
E_new = sum(ci_energy_dressed)
delta_E = dabs(E_new - E_old)
E_past(j) = E_new
j +=1
if(j>4)then
j=1
endif
if(iteration > 4) then
if(delta_E > 1.d-10)then
if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then
print*,'OSCILLATIONS !!!'
oscillations = .True.
i_oscillations +=1
lambda_mrcc_tmp = lambda_mrcc
endif
endif
endif
call save_wavefunction
! if (i_oscillations > 5) then
! exit
! endif
if (iteration > 200) then
exit
endif
print*,'------------'
print*,'VECTOR'
do i = 1, N_det_ref
print*,''
print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
print*,'delta_ii(i,1) = ',delta_ii(i,1)
enddo
print*,'------------'
enddo
call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy")
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
call save_wavefunction
end
subroutine set_generators_bitmasks_as_holes_and_particles
implicit none
integer :: i,k
do k = 1, N_generators_bitmask
do i = 1, N_int
! Pure single part
generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha
generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha
generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta
generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta
! Double excitation
generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha
generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha
generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta
generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta
generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha
generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha
generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta
generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta
enddo
enddo
touch generators_bitmask
end

1
plugins/OVB/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Determinants Psiref_CAS

20
plugins/OVB/README.rst Normal file
View File

@ -0,0 +1,20 @@
=======================
OVB
=======================
The present module proposes an orthogonal Valence Bond analysis
of the wave function, that are the printing of the various Hamiltonian
matrix elements on the basis of the level of ionicity of the components
of the wave function.
Assumptions : it supposes that you have some orthogonal local orbitals within
the active space and that you performed a CI within the active orbitals.
Such CI might be complete or not, no matter.
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

510
plugins/OVB/ovb_components.irp.f Normal file
View File

@ -0,0 +1,510 @@
use bitmasks
BEGIN_PROVIDER [integer, max_number_ionic]
&BEGIN_PROVIDER [integer, min_number_ionic]
BEGIN_DOC
! Maximum and minimum number of ionization in psi_ref
END_DOC
implicit none
integer :: i,j
integer :: n_closed_shell_cas
max_number_ionic = 0
min_number_ionic = 100000
do i = 1, N_det_ref
j = n_closed_shell_cas(psi_ref(1,1,i),n_int)
if(j> max_number_ionic)then
max_number_ionic = j
endif
if(j< min_number_ionic)then
min_number_ionic = j
endif
enddo
print*,'max_number_ionic = ',max_number_ionic
print*,'min_number_ionic = ',min_number_ionic
END_PROVIDER
BEGIN_PROVIDER [integer, ionic_index, (min_number_ionic:max_number_ionic,0:N_det_ref)]
&BEGIN_PROVIDER [double precision, normalization_factor_ionic, (min_number_ionic:max_number_ionic, N_states)]
BEGIN_DOC
! Index of the various determinants in psi_ref according to their level of ionicity
! ionic_index(i,0) = number of determinants in psi_ref having the degree of ionicity "i"
! ionic_index(i,j) = index of the determinants having the degree of ionicity "i"
END_DOC
implicit none
integer :: i,j,k
integer :: n_closed_shell_cas
double precision :: accu
ionic_index = 0
do i = 1, N_det_ref
j = n_closed_shell_cas(psi_ref(1,1,i),n_int)
ionic_index(j,0) +=1
ionic_index(j,ionic_index(j,0)) = i
enddo
do i = min_number_ionic,max_number_ionic
accu = 0.d0
do j = 1, N_states
do k = 1, ionic_index(i,0)
accu += psi_ref_coef_diagonalized(ionic_index(i,k),j) * psi_ref_coef_diagonalized(ionic_index(i,k),j)
enddo
normalization_factor_ionic(i,j) = 1.d0/dsqrt(accu)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, H_OVB_naked, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)]
BEGIN_DOC
! Hamiltonian matrix expressed in the basis of contracted forms in terms of ionic structures
END_DOC
implicit none
integer :: i,j,istate,k,l
double precision :: accu,hij
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
do istate = 1, N_states
accu = 0.d0
do k = 1, ionic_index(i,0)
do l = 1, ionic_index(j,0)
hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
accu += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * &
psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
enddo
enddo
H_OVB_naked(i,j,istate) = accu
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, n_couples_act_orb]
implicit none
n_couples_act_orb = 3
END_PROVIDER
BEGIN_PROVIDER [integer, couples_act_orb, (n_couples_act_orb,2) ]
implicit none
couples_act_orb(1,1) = 20
couples_act_orb(1,2) = 21
couples_act_orb(2,1) = 22
couples_act_orb(2,2) = 23
couples_act_orb(3,1) = 24
couples_act_orb(3,2) = 25
END_PROVIDER
BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom , (n_act_orb,n_act_orb)]
implicit none
BEGIN_DOC
! Hamiltonian matrix elements between the various contracted functions
! that have a negative charge on a given active orbital
END_DOC
integer :: i,j,k,l,jj,ii
integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
double precision :: accu,hij
double precision :: norm
allocate (key_1(N_int,2),key_2(N_int,2))
do i = 1, n_act_orb
j = i ! Diagonal part
norm = 0.d0
accu = 0.d0
do k = 1, n_det_ionic_on_given_atom(i)
norm += psi_coef_mono_ionic_on_given_atom(k,i) **2
do ii = 1, N_int
key_1(ii,1) = psi_det_mono_ionic_on_given_atom(ii,1,k,i)
key_1(ii,2) = psi_det_mono_ionic_on_given_atom(ii,2,k,i)
enddo
do l = 1, n_det_ionic_on_given_atom(j)
do jj = 1, N_int
key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j)
key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j)
enddo
call i_H_j(key_1,key_2,N_int,hij)
accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
enddo
enddo
H_matrix_between_ionic_on_given_atom(i,j) = accu
do j = i+1, n_act_orb ! Extra diagonal part
accu = 0.d0
do k = 1, n_det_ionic_on_given_atom(i)
do jj = 1, N_int
key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i)
key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i)
enddo
do l = 1, n_det_ionic_on_given_atom(j)
do jj = 1, N_int
key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j)
key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j)
enddo
call i_H_j(key_1,key_2,N_int,hij)
accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
enddo
enddo
H_matrix_between_ionic_on_given_atom(i,j) = accu
H_matrix_between_ionic_on_given_atom(j,i) = accu
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom_and_others , (n_act_orb,min_number_ionic:max_number_ionic)]
implicit none
use bitmasks
BEGIN_DOC
! Hamiltonian matrix elements between the various contracted functions
! that have a negative charge on a given active orbital
! and all the other fully contracted OVB structures
END_DOC
integer :: i,j,k,l,jj,ii
integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
double precision :: accu,hij
double precision :: norm
allocate (key_1(N_int,2),key_2(N_int,2))
do i = 1, n_act_orb
do j = min_number_ionic,max_number_ionic
if(j==1)then
H_matrix_between_ionic_on_given_atom_and_others(i,j) = 0.d0
endif
accu = 0.d0
do k = 1, n_det_ionic_on_given_atom(i)
do jj = 1, N_int
key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i)
key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i)
enddo
do l = 1, ionic_index(j,0)
do ii = 1, N_int
key_2(ii,1) = psi_det_ovb(ii,1,l,j)
key_2(ii,2) = psi_det_ovb(ii,2,l,j)
enddo
call i_H_j(key_1,key_2,N_int,hij)
accu += psi_coef_ovb(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
enddo
enddo
H_matrix_between_ionic_on_given_atom_and_others(i,j) = accu
enddo
enddo
print*,'H_matrix_between_ionic_on_given_atom_and_others'
print*,''
do i = 1, n_act_orb
write(*,'(I3,X,100(F16.7))'),H_matrix_between_ionic_on_given_atom_and_others(i,:)
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, n_det_ionic_on_given_atom, (n_act_orb)]
&BEGIN_PROVIDER [double precision, normalization_factor_ionic_on_given_atom, (n_act_orb) ]
&BEGIN_PROVIDER [double precision, psi_coef_mono_ionic_on_given_atom, (N_det_ref,n_act_orb) ]
&BEGIN_PROVIDER [integer(bit_kind), psi_det_mono_ionic_on_given_atom, (N_int,2,N_det_ref,n_act_orb)]
implicit none
use bitmasks
BEGIN_DOC
! number of determinants that are mono ionic with the negative charge
! on a given atom, normalization_factor, array of determinants,and coefficients
END_DOC
integer :: i,j,k,l
ionicity_level = 1
integer :: ionicity_level
logical :: doubly_occupied_array(n_act_orb)
n_det_ionic_on_given_atom = 0
normalization_factor_ionic_on_given_atom = 0.d0
do i = 1, ionic_index(ionicity_level,0)
call give_index_of_doubly_occ_in_active_space(psi_det(1,1,ionic_index(ionicity_level,i)),doubly_occupied_array)
do j = 1, n_act_orb
if(doubly_occupied_array(j))then
n_det_ionic_on_given_atom(j) += 1
normalization_factor_ionic_on_given_atom(j) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2
do k = 1, N_int
psi_det_mono_ionic_on_given_atom(k,1,n_det_ionic_on_given_atom(j),j) = psi_det(k,1,ionic_index(ionicity_level,i))
psi_det_mono_ionic_on_given_atom(k,2,n_det_ionic_on_given_atom(j),j) = psi_det(k,2,ionic_index(ionicity_level,i))
enddo
psi_coef_mono_ionic_on_given_atom(n_det_ionic_on_given_atom(j),j) = psi_ref_coef_diagonalized(ionic_index(1,i),1)
endif
enddo
enddo
integer :: i_count
i_count = 0
do j = 1, n_act_orb
i_count += n_det_ionic_on_given_atom(j)
normalization_factor_ionic_on_given_atom(j) = 1.d0/dsqrt(normalization_factor_ionic_on_given_atom(j))
enddo
if(i_count.ne.ionic_index(ionicity_level,0))then
print*,'PB with n_det_ionic_on_given_atom'
print*,'i_count = ',i_count
print*,'ionic_index(ionicity_level,0)',ionic_index(ionicity_level,0)
stop
endif
do j = 1, n_act_orb
do i = 1, n_det_ionic_on_given_atom(j)
psi_coef_mono_ionic_on_given_atom(i,j) = psi_coef_mono_ionic_on_given_atom(i,j) * normalization_factor_ionic_on_given_atom(j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [integer(bit_kind), psi_det_ovb, (N_int,2,N_det_ref,min_number_ionic:max_number_ionic)]
&BEGIN_PROVIDER [double precision, psi_coef_ovb, (N_det_ref,min_number_ionic:max_number_ionic) ]
implicit none
BEGIN_DOC
! Array of the determinants belonging to each ovb structures (neutral, mono ionic, bi ionic etc ...)
! together with the arrays of coefficients
END_DOC
integer :: i,j,k,l
use bitmasks
integer :: ionicity_level,i_count
double precision :: accu
do ionicity_level = min_number_ionic,max_number_ionic
accu = 0.d0
do i = 1, ionic_index(ionicity_level,0)
do j = 1, N_int
psi_det_ovb(j,1,i,ionicity_level) = psi_det(j,1,ionic_index(ionicity_level,i))
psi_det_ovb(j,2,i,ionicity_level) = psi_det(j,2,ionic_index(ionicity_level,i))
enddo
psi_coef_ovb(i,ionicity_level) = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) * normalization_factor_ionic(ionicity_level,1)
accu += psi_coef_ovb(i,ionicity_level)**2
enddo
accu = 1.d0/dsqrt(accu)
do i = 1, ionic_index(ionicity_level,0)
psi_coef_ovb(i,ionicity_level) = psi_coef_ovb(i,ionicity_level) * accu
enddo
accu = 0.d0
do i = 1, ionic_index(ionicity_level,0)
accu += psi_coef_ovb(i,ionicity_level) **2
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, H_matrix_psi_det_ovb, (min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)]
implicit none
BEGIN_DOC
! H matrix between the fully contracted OVB forms
END_DOC
integer :: i,j,k,l,jj,ii
integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
use bitmasks
double precision :: accu,hij
double precision :: norm
allocate (key_1(N_int,2),key_2(N_int,2))
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
accu = 0.d0
do k = 1, ionic_index(i,0)
do ii = 1, N_int
key_1(ii,1) = psi_det_ovb(ii,1,k,i)
key_1(ii,2) = psi_det_ovb(ii,2,k,i)
enddo
do l = 1, ionic_index(j,0)
do ii = 1, N_int
key_2(ii,1) = psi_det_ovb(ii,1,l,j)
key_2(ii,2) = psi_det_ovb(ii,2,l,j)
enddo
call i_H_j(key_1,key_2,N_int,hij)
accu += psi_coef_ovb(l,j) * psi_coef_ovb(k,i) * hij
enddo
enddo
H_matrix_psi_det_ovb(i,j) = accu
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, number_first_ionic_couples]
&BEGIN_PROVIDER [logical , is_a_first_ionic_couple, (N_det_ref)]
&BEGIN_PROVIDER [double precision, normalization_factor_special_first_ionic, (2)]
implicit none
BEGIN_DOC
! Number of determinants belonging to the class of first ionic
! AND that have a couple of positive/negative charge belonging
! to a couple of orbital couples_act_orb
! If is_a_first_ionic_couple(i) = .True. then this determinant is a first ionic
! and have a couple of positive/negative charge belonging
! to a couple of orbital couples_act_orb
! normalization factor (1) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .True.)
! normalization factor (2) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .False.)
END_DOC
integer :: i,j
use bitmasks
number_first_ionic_couples = 0
integer :: ionicity_level
logical :: couples_out(0:n_couples_act_orb)
integer(bit_kind) :: key_tmp(N_int,2)
ionicity_level = 1
normalization_factor_special_first_ionic = 0.d0
do i = 1, ionic_index(ionicity_level,0)
do j = 1, N_int
key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i))
key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i))
enddo
call doubly_occ_empty_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out)
if(couples_out(0))then
number_first_ionic_couples +=1
is_a_first_ionic_couple(i) = .True.
normalization_factor_special_first_ionic(1) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2
else
is_a_first_ionic_couple(i) = .False.
normalization_factor_special_first_ionic(2) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2
endif
enddo
normalization_factor_special_first_ionic(1) = 1.d0/dsqrt(normalization_factor_special_first_ionic(1))
normalization_factor_special_first_ionic(2) = 1.d0/dsqrt(normalization_factor_special_first_ionic(2))
print*,'number_first_ionic_couples = ',number_first_ionic_couples
END_PROVIDER
BEGIN_PROVIDER [integer, number_neutral_no_hund_couples]
&BEGIN_PROVIDER [logical , is_a_neutral_no_hund_couple, (N_det_ref)]
&BEGIN_PROVIDER [double precision, normalization_factor_neutra_no_hund_couple, (2)]
&BEGIN_PROVIDER [double precision, ratio_hund_no_hund ]
implicit none
BEGIN_DOC
! Number of determinants belonging to the class of neutral determinants
! AND that have a couple of alpha beta electrons in couple of orbital couples_act_orb
! If is_a_neutral_no_hund_couple(i) = .True. then this determinant is a neutral determinants
! and have a a couple of alpha beta electrons in couple of orbital couples_act_orb
! normalization factor (1) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .True.)
! normalization factor (2) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .False.)
END_DOC
integer :: i,j
use bitmasks
number_neutral_no_hund_couples = 0
integer :: ionicity_level
logical :: couples_out(0:n_couples_act_orb)
integer(bit_kind) :: key_tmp(N_int,2)
integer :: ifirst_hund,ifirst_no_hund
double precision :: coef_ref_hund,coef_ref_no_hund
ifirst_hund = 0
ifirst_no_hund = 0
ionicity_level = 0
normalization_factor_neutra_no_hund_couple = 0.d0
do i = 1, ionic_index(ionicity_level,0)
do j = 1, N_int
key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i))
key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i))
enddo
call neutral_no_hund_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out)
if(couples_out(0))then
if(ifirst_no_hund == 0)then
coef_ref_no_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1)
ifirst_no_hund = 1
endif
number_neutral_no_hund_couples +=1
is_a_neutral_no_hund_couple(i) = .True.
normalization_factor_neutra_no_hund_couple(1) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2
else
if(ifirst_hund == 0)then
coef_ref_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1)
ifirst_hund = 1
endif
is_a_neutral_no_hund_couple(i) = .False.
normalization_factor_neutra_no_hund_couple(2) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2
endif
enddo
ratio_hund_no_hund = coef_ref_no_hund/coef_ref_hund
normalization_factor_neutra_no_hund_couple(1) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(1))
normalization_factor_neutra_no_hund_couple(2) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(2))
print*,'number_neutral_no_hund_couples = ',number_neutral_no_hund_couples
END_PROVIDER
BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic, (2,min_number_ionic:max_number_ionic,n_states)]
&BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic_between_ionic, (2,2,n_states)]
BEGIN_DOC
! H_OVB_naked_first_ionic(1,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
! and the contracted ith ionic form
! if i == 1 not defined
! H_OVB_naked_first_ionic(2,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False
! and the contracted ith ionic form
! if i == 1 not defined
! H_OVB_naked_first_ionic_between_ionic(1,1) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
! and the first ionic determinants belonging to is_a_first_ionic_couple = True
! H_OVB_naked_first_ionic_between_ionic(1,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
! and the first ionic determinants belonging to is_a_first_ionic_couple = False
! H_OVB_naked_first_ionic_between_ionic(2,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False
! and the first ionic determinants belonging to is_a_first_ionic_couple = False
END_DOC
implicit none
integer :: i,j,istate,k,l
double precision :: accu_1,accu_2,hij
H_OVB_naked_first_ionic = 0.d0
H_OVB_naked_first_ionic_between_ionic = 0.d0
i = 1
do j = min_number_ionic,max_number_ionic
if(j==1)cycle
do istate = 1, N_states
accu_1 = 0.d0
accu_2 = 0.d0
do k = 1, ionic_index(i,0)
if(is_a_first_ionic_couple(k))then
do l = 1, ionic_index(j,0)
hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
accu_1 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(1) * &
psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
enddo
else
do l = 1, ionic_index(j,0)
hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
accu_2 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(2) * &
psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
enddo
endif
enddo
H_OVB_naked_first_ionic(1,j,istate) = accu_1
H_OVB_naked_first_ionic(2,j,istate) = accu_2
enddo
enddo
do istate = 1, N_states
accu_1 = 0.d0
accu_2 = 0.d0
integer :: i_count
i_count = 0
do k = 1, ionic_index(1,0)
do l = 1, ionic_index(1,0)
hij = ref_hamiltonian_matrix(ionic_index(1,k),ionic_index(1,l))
accu_1 = hij * psi_ref_coef_diagonalized(ionic_index(1,k),istate) * psi_ref_coef_diagonalized(ionic_index(1,l),istate)
if(is_a_first_ionic_couple(k).and. is_a_first_ionic_couple(l))then
H_OVB_naked_first_ionic_between_ionic(1,1,istate) += accu_1 * normalization_factor_special_first_ionic(1) **2
elseif(is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then
i_count += 1
H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * &
normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2)
! elseif(is_a_first_ionic_couple(l).and. .not.is_a_first_ionic_couple(k))then
! i_count += 1
! H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * &
! normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2)
elseif(.not.is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then
H_OVB_naked_first_ionic_between_ionic(2,2,istate) += accu_1 * normalization_factor_special_first_ionic(2) **2
endif
enddo
enddo
enddo
print*,'i_count = ',i_count
print*,'number_first_ionic_couples**2 = ',ionic_index(1,0) * number_first_ionic_couples
double precision :: convert_hartree_ev
convert_hartree_ev = 27.211399d0
print*,'Special H matrix'
do i = 1,2
write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic(i,:,1)
enddo
print*,'Special H matrix bis'
do i = 1,2
write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic_between_ionic(i,:,1)
enddo
END_PROVIDER

27
plugins/OVB/print_ovb.irp.f Normal file
View File

@ -0,0 +1,27 @@
program print_OVB
implicit none
read_wf = .True.
call provide_all
end
subroutine provide_all
implicit none
integer :: i,j,k,l,istate
do istate= 1, N_states
print*,'-------------------'
print*,'ISTATE = ',istate
print*,'-------------------'
print*,'CAS MATRIX '
print*,''
do i = min_number_ionic,max_number_ionic
write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_naked(i,:,istate)
enddo
print*,''
print*,'-------------------'
print*,'-------------------'
enddo
end

1
plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Determinants

345
plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f Normal file
View File

@ -0,0 +1,345 @@
use bitmasks
BEGIN_PROVIDER [integer, mo_inp_num]
implicit none
BEGIN_DOC
! This is the number of orbitals involved in the entanglement calculation.
! It is taken equal to the number of active orbitals n_act_orb.
END_DOC
mo_inp_num = n_act_orb
END_PROVIDER
BEGIN_PROVIDER [integer, mo_inp_list, (N_int*bit_kind_size)]
&BEGIN_PROVIDER [integer, mo_inp_list_rev, (mo_tot_num)]
&BEGIN_PROVIDER [integer(bit_kind), mo_inp_bit_list, (N_int)]
implicit none
BEGIN_DOC
! mo_inp_list is the list of the orbitals involved in the entanglement calculation.
! It is taken equal to the list of active orbitals list_act.
! mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
END_DOC
integer :: i
do i = 1, mo_inp_num
mo_inp_list(i)=list_act(i)
enddo
do i = 1, mo_inp_num
mo_inp_list_rev(mo_inp_list(i))=i
enddo
call list_to_bitstring( mo_inp_bit_list, mo_inp_list, mo_inp_num, N_int)
END_PROVIDER
BEGIN_PROVIDER [double precision, entropy_one_orb, (mo_inp_num)]
&BEGIN_PROVIDER [double precision, entropy_two_orb, (mo_inp_num,mo_inp_num)]
implicit none
BEGIN_DOC
! entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
! entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
END_DOC
double precision, allocatable :: ro1(:,:),ro2(:,:,:)
integer :: i,j,k,l,ii,jj,iii,istate,kl,info
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha, n_occ_beta
logical, allocatable :: zocc(:,:)
logical :: zalpha, zbeta, zalpha2, zbeta2
integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,spin1,spin2
double precision :: phase
integer(bit_kind) :: key_tmp(N_int), key_tmp2(N_int)
integer :: ip
double precision, parameter :: eps=10.d0**(-14)
double precision :: w(16), work(3*16)
allocate(ro1(4,mo_inp_num),ro2(16,16,(mo_inp_num*(mo_inp_num-1)/2)))
entropy_one_orb = 0.d0
entropy_two_orb = 0.d0
ro1 = 0.d0
ro2 = 0.d0
allocate (occ(N_int*bit_kind_size,2))
allocate (zocc(mo_tot_num,2))
istate = 1 !Only GS, to be generalized...
do ii=1,N_det
! We get the occupation of the alpha electrons in occ(:,1)
call bitstring_to_list(psi_det(1,1,ii), occ(1,1), n_occ_alpha, N_int)
! We get the occupation of the beta electrons in occ(:,2)
call bitstring_to_list(psi_det(1,2,ii), occ(1,2), n_occ_beta, N_int)
zocc = .false.
do i=1,n_occ_alpha
zocc(occ(i,1),1)=.true.
enddo
do i=1,n_occ_beta
zocc(occ(i,2),2)=.true.
enddo
do k=1,mo_inp_num
zalpha = zocc(mo_inp_list(k),1)
zbeta = zocc(mo_inp_list(k),2)
! mono start
if (zbeta.and.zalpha) then
ro1(4,k) = ro1(4,k) + psi_coef(ii,istate)**2 ! double occupied
elseif (zalpha) then
ro1(2,k) = ro1(2,k) + psi_coef(ii,istate)**2 ! single alpha
elseif (zbeta) then
ro1(3,k) = ro1(3,k) + psi_coef(ii,istate)**2 ! single beta
else
ro1(1,k) = ro1(1,k) + psi_coef(ii,istate)**2 ! empty
endif
! mono stop
! double start
if (k.eq.mo_inp_num) cycle
do l=k+1,mo_inp_num
kl=(l-1)*(l-2)/2+k
zalpha2 = zocc(mo_inp_list(l),1)
zbeta2 = zocc(mo_inp_list(l),2)
if (zbeta.and.zalpha.and.zbeta2.and.zalpha2) then
ro2(16,16,kl) = ro2(16,16,kl) + psi_coef(ii,istate)**2 ! both double occupied
else if (zbeta.and.zalpha.and.zbeta2) then
ro2(15,15,kl) = ro2(15,15,kl) + psi_coef(ii,istate)**2 ! one double, one beta
else if (zbeta.and.zalpha.and.zalpha2) then
ro2(13,13,kl) = ro2(13,13,kl) + psi_coef(ii,istate)**2 ! one double, one alpha
else if (zbeta.and.zbeta2.and.zalpha2) then
ro2(14,14,kl) = ro2(14,14,kl) + psi_coef(ii,istate)**2 ! one beta, one double
else if (zalpha.and.zbeta2.and.zalpha2) then
ro2(12,12,kl) = ro2(12,12,kl) + psi_coef(ii,istate)**2 ! one alpha, one double
else if (zalpha.and.zbeta) then
ro2(11,11,kl) = ro2(11,11,kl) + psi_coef(ii,istate)**2 ! one double, one empty
else if (zbeta2.and.zalpha2) then
ro2(8,8,kl) = ro2(8,8,kl) + psi_coef(ii,istate)**2 ! one empty, one double
else if (zbeta.and.zalpha2) then
ro2(10,10,kl) = ro2(10,10,kl) + psi_coef(ii,istate)**2 ! one beta, one alpha
else if (zalpha.and.zbeta2) then
ro2(9,9,kl) = ro2(9,9,kl) + psi_coef(ii,istate)**2 ! one alpha, one beta
else if (zbeta.and.zbeta2) then
ro2(7,7,kl) = ro2(7,7,kl) + psi_coef(ii,istate)**2 ! one beta, one beta
else if (zalpha.and.zalpha2) then
ro2(6,6,kl) = ro2(6,6,kl) + psi_coef(ii,istate)**2 ! one alpha, one alpha
else if (zbeta) then
ro2(5,5,kl) = ro2(5,5,kl) + psi_coef(ii,istate)**2 ! one beta, one empty
else if (zbeta2) then
ro2(4,4,kl) = ro2(4,4,kl) + psi_coef(ii,istate)**2 ! one empty, one beta
else if (zalpha) then
ro2(3,3,kl) = ro2(3,3,kl) + psi_coef(ii,istate)**2 ! one alpha, one empty
else if (zalpha2) then
ro2(2,2,kl) = ro2(2,2,kl) + psi_coef(ii,istate)**2 ! one empty, one alpha
else
ro2(1,1,kl) = ro2(1,1,kl) + psi_coef(ii,istate)**2 ! both empty
end if
enddo
enddo
! stop double
if (ii.eq.N_det) cycle
!Off Diagonal Elements
do jj=ii+1,N_det
call get_excitation_degree(psi_det(1,1,ii),psi_det(1,1,jj),degree,N_int)
if (degree.gt.2) cycle
ip=0
do iii =1,N_int
key_tmp(iii) = ior(xor(psi_det(iii,1,ii),psi_det(iii,1,jj)),xor(psi_det(iii,2,ii),psi_det(iii,2,jj)))
ip += popcnt(key_tmp(iii))
enddo
if (ip.ne.2) cycle !They involve more than 2 orbitals.
ip=0
do iii=1,N_int
ip += popcnt(iand(key_tmp(iii),mo_inp_bit_list(iii)))
enddo
if (ip.ne.2) cycle !They do not involve orbitals of the list.
if (degree.eq.2) then
call get_double_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2)
k=mo_inp_list_rev(h1)
l=mo_inp_list_rev(p1)
if (k.gt.l) then
kl=(k-1)*(k-2)/2+l
else
kl=(l-1)*(l-2)/2+k
endif
if ((.not.zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))&
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
ro2(8,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(11,8,kl) = ro2(8,11,kl)
endif
if ((zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))&
.and.(.not.zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(10,9,kl) = ro2(9,10,kl)
endif
if ((zocc(mo_inp_list(k),1)).and.(.not.zocc(mo_inp_list(k),2))&
.and.(.not.zocc(mo_inp_list(l),1)).and.(zocc(mo_inp_list(l),2))) then
ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(10,9,kl) = ro2(9,10,kl)
endif
endif
if (degree.eq.1) then
call get_mono_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2)
k=mo_inp_list_rev(h1)
l=mo_inp_list_rev(p1)
if (k.gt.l) then
kl=(k-1)*(k-2)/2+l
else
kl=(l-1)*(l-2)/2+k
endif
if ((.not.(zocc(mo_inp_list(l),2))).and.&
(.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),1))&
.and.(.not.(zocc(mo_inp_list(k),2)))) then
ro2(2,3,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(3,2,kl)=ro2(2,3,kl)
endif
if ((.not.(zocc(mo_inp_list(l),2))).and.&
(.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),2))&
.and.(.not.(zocc(mo_inp_list(k),1)))) then
ro2(4,5,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(5,4,kl)=ro2(4,5,kl)
endif
if ((.not.(zocc(mo_inp_list(l),1))).and.& !k doubly occupied, l empty
(.not.(zocc(mo_inp_list(l),2))).and.&
(zocc(mo_inp_list(k),1)).and.&
(zocc(mo_inp_list(k),2))) then
if (k.gt.l) then
if (spin1.eq.1) then !spin alpha
ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(9,8,kl)=ro2(8,9,kl)
else !spin beta
ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(10,8,kl)=ro2(8,10,kl)
endif
else ! k.lt.l
if (spin1.eq.1) then !spin alpha
ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(11,10,kl)=ro2(10,11,kl)
else !spin beta
ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(11,9,kl)=ro2(9,11,kl)
endif
endif
endif
if ((.not.(zocc(mo_inp_list(l),1))).and.& !k alpha, l beta
(.not.(zocc(mo_inp_list(k),2))).and.&
(zocc(mo_inp_list(k),1)).and.&
(zocc(mo_inp_list(l),2))) then
if (k.gt.l) then
if (spin1.eq.1) then !spin alpha
ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(11,10,kl)=ro2(10,11,kl)
else !spin beta
print*, "problem in k alpha l beta k.gt.l spin beta"
endif
else ! k.lt.l
if (spin1.eq.1) then !spin alpha
ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(9,8,kl)=ro2(8,9,kl)
else !spin beta
print*, "problem in k alpha l beta k.lt.l spin beta"
endif
endif
endif
if ((.not.(zocc(mo_inp_list(k),1))).and.& !k beta, l alpha
(.not.(zocc(mo_inp_list(l),2))).and.&
(zocc(mo_inp_list(l),1)).and.&
(zocc(mo_inp_list(k),2))) then
if (k.gt.l) then
if (spin1.eq.2) then !spin beta
ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(11,9,kl)=ro2(9,11,kl)
else !spin alpha
print*, "problem in k beta l alpha k.gt.l spin alpha"
endif
else ! k.lt.l
if (spin1.eq.2) then !spin beta
ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
ro2(10,8,kl)=ro2(8,10,kl)
else !spin alpha
print*, "problem in k beta l alpha k.lt.l spin alpha"
endif
endif
endif
if (zocc(mo_inp_list(l),1).and.(.not.zocc(mo_inp_list(l),2))&
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
ro2(12,13,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(13,12,kl) = ro2(12,13,kl)
endif
if (zocc(mo_inp_list(l),2).and.(.not.zocc(mo_inp_list(l),1))&
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
ro2(14,15,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate)
ro2(15,14,kl) = ro2(14,15,kl)
endif
endif
enddo
enddo
entropy_one_orb=0.d0
do k=1,mo_inp_num
do i=1,4
if (ro1(i,k).ge.eps) then
entropy_one_orb(k) = entropy_one_orb(k)-ro1(i,k)*log(ro1(i,k))
endif
enddo
enddo
entropy_two_orb=0.d0
do k=1,mo_inp_num
do l=1,k
if (k.eq.l) cycle
kl=(k-1)*(k-2)/2+l
call dsyev('N','U',16,ro2(1,1,kl),16,w,work,3*16,info)
if (info.ne.0) then
write(*,*) "Errore in dsyev"
endif
do j=1,16
if (w(j).ge.eps) then
entropy_two_orb(k,l) = entropy_two_orb(k,l)-w(j)*log(w(j))
entropy_two_orb(l,k) = entropy_two_orb(k,l)
elseif ((w(j)).lt.(-eps)) then
write(6,*) "Negative Eigenvalue. You have a big problem..."
write(6,*) w(j)
stop
endif
enddo
enddo
enddo
deallocate (occ,zocc,ro1,ro2)
END_PROVIDER
BEGIN_PROVIDER [double precision, mutinf, (mo_inp_num,mo_inp_num)]
implicit none
BEGIN_DOC
!mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij]
!see the refence: 10.1016/j.chemphys.2005.10.018
END_DOC
integer :: i,j
! mutal information:
mutinf = 0.d0
do i=1,mo_inp_num
do j=1,mo_inp_num
if (j.eq.i) cycle
mutinf(i,j)=-0.5d0*(entropy_two_orb(i,j)-entropy_one_orb(i)-entropy_one_orb(j))
enddo
enddo
END_PROVIDER

65
plugins/Orbital_Entanglement/README.rst Normal file
View File

@ -0,0 +1,65 @@
====================
Orbital_Entanglement
====================
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
.. image:: tree_dependency.png
* `Determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants>`_
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
`entropy_one_orb <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L35>`_
entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
`entropy_two_orb <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L36>`_
entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
`mo_inp_bit_list <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L15>`_
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
It is taken equal to the list of active orbitals list_act.
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
`mo_inp_list <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L13>`_
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
It is taken equal to the list of active orbitals list_act.
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
`mo_inp_list_rev <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L14>`_
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
It is taken equal to the list of active orbitals list_act.
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
`mo_inp_num <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L3>`_
This is the number of orbitals involved in the entanglement calculation.
It is taken equal to the number of active orbitals n_act_orb.
`mutinf <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L330>`_
mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij]
see the refence: 10.1016/j.chemphys.2005.10.018
`pouet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/print_entanglement.irp.f#L1>`_
Undocumented
`routine <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/print_entanglement.irp.f#L9>`_
Undocumented

46
plugins/Orbital_Entanglement/print_entanglement.irp.f Normal file
View File

@ -0,0 +1,46 @@
program pouet
implicit none
read_wf = .true.
touch read_wf
call routine
end
subroutine routine
implicit none
integer:: i,j
write(6,*) 'Total orbitals: ', mo_tot_num
write(6,*) 'Orbitals for entanglement calculation: ', mo_inp_num
write(6,*) 'Index: ',(mo_inp_list(i),i=1,mo_inp_num)
write(6,*)
write(6,*) "s1: One-Orbital von Neumann entropy"
write(6,'(1000f8.5)') entropy_one_orb
write(6,*)
write(6,*) "s2: Two-Orbital von Neumann entropy"
do i=1,mo_inp_num
write(6,'(1000f8.5)') (entropy_two_orb(i,j),j=1,mo_inp_num)
enddo
write(6,*)
! mutal information:
write(6,*) "Mutual Information (Entanglement)"
do i=1,mo_inp_num
write(6,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num)
enddo
open(17,file=(trim(ezfio_filename)//".entanglement"),status='unknown',form='formatted')
write(17,'(1000f8.5)') entropy_one_orb
do i=1,mo_inp_num
write(17,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num)
enddo
close(17)
write(6,*)
write(6,*) "The .entanglement file is the input for the entanglement.py script."
write(6,*) "You can find the script in the directory Scripts of QP."
write(6,*)
end

21
plugins/Orbital_Entanglement/tree_dependency Normal file
View File

@ -0,0 +1,21 @@
// ['Orbital_Entanglement']
digraph {
Orbital_Entanglement [fontcolor=red]
Orbital_Entanglement -> Determinants
Determinants -> Integrals_Monoelec
Integrals_Monoelec -> MO_Basis
MO_Basis -> AO_Basis
AO_Basis -> Nuclei
Nuclei -> Ezfio_files
Nuclei -> Utils
MO_Basis -> Electrons
Electrons -> Ezfio_files
Integrals_Monoelec -> Pseudo
Pseudo -> Nuclei
Determinants -> Integrals_Bielec
Integrals_Bielec -> Pseudo
Integrals_Bielec -> Bitmask
Bitmask -> MO_Basis
Integrals_Bielec -> ZMQ
ZMQ -> Utils
}

0
plugins/Orbital_Entanglement/tree_dependency.png Normal file
View File

135
plugins/Properties/hyperfine_constants.irp.f Normal file
View File

@ -0,0 +1,135 @@
BEGIN_PROVIDER [double precision, spin_density_at_nucleous, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
do k = 1, ao_num
accu += one_body_spin_density_ao(k,j) * aos_array(k) * aos_array(j)
enddo
enddo
spin_density_at_nucleous(i) = accu
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_from_mo, (nucl_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_per_mo, (nucl_num,mo_tot_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: contrib
double precision :: mo_values(mo_tot_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
spin_density_at_nucleous_from_mo(i) = 0.d0
do k = 1, mo_tot_num
mo_values(k) = 0.d0
do j = 1, ao_num
mo_values(k) += mo_coef(j,k) * aos_array(j)
enddo
enddo
do k = 1, mo_tot_num
spin_density_at_nucleous_contrib_per_mo(i,k) = 0.d0
do m = 1, mo_tot_num
contrib = one_body_spin_density_mo(k,m) * mo_values(k) * mo_values(m)
spin_density_at_nucleous_from_mo(i) += contrib
spin_density_at_nucleous_contrib_per_mo(i,k) += contrib
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo, (nucl_num,mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo_test, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
spin_density_at_nucleous_contrib_mo_test = 0.d0
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: c_i1,c_j1
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do k = 1, mo_tot_num
do m = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
c_i1 = mo_coef(j,k)
do l = 1, ao_num
c_j1 = c_i1*mo_coef(l,m)
accu += one_body_spin_density_mo(k,m) * aos_array(l) * aos_array(j) * c_j1
enddo
enddo
spin_density_at_nucleous_contrib_mo(i,k,m) = accu
spin_density_at_nucleous_contrib_mo_test(i) += accu
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, conversion_factor_mhz_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_gauss_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_cm_1_hcc, (100)]
BEGIN_DOC
! Conversion factor for the calculation of the hcc, according to the nuclear charge
END_DOC
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
! hydrogen
conversion_factor_mhz_hcc(1) = 4469.84692227102460d0
conversion_factor_gauss_hcc(1) = 1594.95296390862904d0
conversion_factor_cm_1_hcc(1) = 1490.98044430157870d0
! Li
conversion_factor_mhz_hcc(3) = 1737.2746512855997d0
conversion_factor_gauss_hcc(3) = 619.9027742370165d0
conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
! carbon
conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
conversion_factor_gauss_hcc(6) = 401.136570647523058d0
conversion_factor_cm_1_hcc(6) = 374.987097339830086d0
! nitrogen
conversion_factor_mhz_hcc(7) = 323.102093833793390d0
conversion_factor_gauss_hcc(7) = 115.290892768082614d0
conversion_factor_cm_1_hcc(7) = 107.775257586297698d0
! Oxygen
conversion_factor_mhz_hcc(8) = -606.1958551736545d0
conversion_factor_gauss_hcc(8) = -216.30574771560407d0
conversion_factor_cm_1_hcc(8) = -202.20517197179822d0
END_PROVIDER
BEGIN_PROVIDER [double precision, iso_hcc_mhz, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_gauss, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_cm_1, (nucl_num)]
BEGIN_DOC
! isotropic hyperfine coupling constants among the various atoms
END_DOC
integer :: i
do i = 1, nucl_num
iso_hcc_mhz(i) = conversion_factor_mhz_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !* 0.5d0
iso_hcc_gauss(i) = conversion_factor_gauss_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)!* 0.5d0
iso_hcc_cm_1(i) = conversion_factor_cm_1_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !*0.5d0
enddo
END_PROVIDER

107
plugins/Properties/mulliken.irp.f Normal file
View File

@ -0,0 +1,107 @@
BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
spin_population = 0.d0
do i = 1, ao_num
do j = 1, ao_num
spin_population(j,i) = one_body_spin_density_ao(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)]
implicit none
integer :: i
double precision :: accu
spin_population_angular_momentum = 0.d0
do i = 1, ao_num
spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_gross_orbital_product, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product for the spin population
END_DOC
do i = 1, ao_num
do j = 1, ao_num
spin_gross_orbital_product(i) += spin_population(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_spin_densities, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!ATOMIC SPIN POPULATION (ALPHA MINUS BETA)
END_DOC
mulliken_spin_densities = 0.d0
do i = 1, ao_num
mulliken_spin_densities(ao_nucl(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, electronic_population_alpha, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [double precision, electronic_population_beta, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
electronic_population_alpha = 0.d0
electronic_population_beta = 0.d0
do i = 1, ao_num
do j = 1, ao_num
electronic_population_alpha(j,i) = one_body_dm_ao_alpha(i,j) * ao_overlap(i,j)
electronic_population_beta(j,i) = one_body_dm_ao_beta(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, gross_orbital_product_alpha, (ao_num)]
&BEGIN_PROVIDER [double precision, gross_orbital_product_beta, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product
END_DOC
do i = 1, ao_num
do j = 1, ao_num
gross_orbital_product_alpha(i) += electronic_population_alpha(j,i)
gross_orbital_product_beta(i) += electronic_population_beta(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_densities_alpha, (nucl_num)]
&BEGIN_PROVIDER [double precision, mulliken_densities_beta, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!
END_DOC
mulliken_densities_alpha = 0.d0
mulliken_densities_beta = 0.d0
do i = 1, ao_num
mulliken_densities_alpha(ao_nucl(i)) += gross_orbital_product_alpha(i)
mulliken_densities_beta(ao_nucl(i)) += gross_orbital_product_beta(i)
enddo
END_PROVIDER

17
plugins/Properties/print_hcc.irp.f Normal file
View File

@ -0,0 +1,17 @@
program print_hcc
implicit none
read_wf = .True.
touch read_wf
call test
end
subroutine test
implicit none
double precision :: accu
integer :: i,j
print*,'Z AU GAUSS MHZ cm^-1'
do i = 1, nucl_num
write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
enddo
end

35
plugins/Properties/print_mulliken.irp.f Normal file
View File

@ -0,0 +1,35 @@
program print_mulliken
implicit none
read_wf = .True.
touch read_wf
print*,'Mulliken spin densities'
call test
end
subroutine test
double precision :: accu
integer :: i
integer :: j
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += spin_gross_orbital_product(i)
write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += spin_population_angular_momentum(i)
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu
enddo
end

21
plugins/Psiref_Utils/psi_ref_utils.irp.f
View File

@ -125,7 +125,7 @@ BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)]
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)]
BEGIN_PROVIDER [double precision, psi_ref_coef_diagonalized, (N_det_ref,N_states)]
&BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)]
implicit none
integer :: i,j
@ -137,9 +137,11 @@ END_PROVIDER
do i = 1, N_states
psi_ref_energy_diagonalized(i) = eigenvalues(i)
do j = 1, N_det_ref
psi_coef_ref_diagonalized(j,i) = eigenvectors(j,i)
psi_ref_coef_diagonalized(j,i) = eigenvectors(j,i)
enddo
enddo
deallocate (eigenvectors)
deallocate (eigenvalues)
END_PROVIDER
@ -264,3 +266,18 @@ integer function get_index_in_psi_ref_sorted_bit(key,Nint)
end
BEGIN_PROVIDER [double precision, ref_hamiltonian_matrix, (n_det_ref,n_det_ref)]
BEGIN_DOC
! H matrix in the Reference space
END_DOC
implicit none
integer :: i,j
double precision :: hij
do i = 1, N_det_ref
do j = 1, N_det_ref
call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
ref_hamiltonian_matrix(i,j) = hij
enddo
enddo
END_PROVIDER

232
plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py
View File

@ -1,232 +0,0 @@
#!/usr/bin/python
print "#QP -> QMCPACK"
from ezfio import ezfio
import sys
ezfio_path = sys.argv[1]
ezfio.set_file(ezfio_path)
do_pseudo = ezfio.get_pseudo_do_pseudo()
if do_pseudo:
print "do_pseudo True"
zcore = ezfio.get_pseudo_nucl_charge_remove()
else:
print "do_pseudo False"
n_det =ezfio.get_determinants_n_det()
if n_det == 1:
print "multi_det False"
else:
print "multi_det True"
ao_num = ezfio.get_ao_basis_ao_num()
print "ao_num", ao_num
mo_num = ezfio.get_mo_basis_mo_tot_num()
print "mo_num", mo_num
alpha = ezfio.get_electrons_elec_alpha_num()
beta = ezfio.get_electrons_elec_beta_num()
print "elec_alpha_num", alpha
print "elec_beta_num", beta
print "elec_tot_num", alpha + beta
print "spin_multiplicity", 2*(alpha-beta)+1
l_label = ezfio.get_nuclei_nucl_label()
l_charge = ezfio.get_nuclei_nucl_charge()
l_coord = ezfio.get_nuclei_nucl_coord()
l_coord_str = [" ".join(map(str,i)) for i in l_coord]
print "nucl_num",len(l_label)
print "Atomic coord in Bohr"
for i,t in enumerate(zip(l_label,l_charge,l_coord_str)):
try :
l = (t[0],t[1]+zcore[i],t[1])
except NameError:
l = t
print " ".join(map(str,l))
import subprocess
process = subprocess.Popen(['qp_print_basis', ezfio_path], stdout=subprocess.PIPE)
out, err = process.communicate()
basis_raw, sym_raw, mo_raw = out.split("\n\n\n")
basis_without_header = "\n".join(basis_raw.split("\n")[7:])
for i,l in enumerate(l_label):
basis_without_header=basis_without_header.replace('Atom {0}'.format(i+1),l)
print "BEGIN_BASIS_SET"
print ""
print basis_without_header
print "END_BASIS_SET"
# _
# |\/| / \ _
# | | \_/ _>
#
def same_character(item1):
return item1==item1[0]* len(item1)
def compare_gamess_style(item1, item2):
if len(item1) < len(item2):
return -1
elif len(item1) > len(item2):
return 1
elif same_character(item1) and same_character(item2):
if item1 < item2:
return -1
else:
return 1
elif same_character(item1) and not same_character(item2):
return -1
elif not same_character(item1) and same_character(item2):
return 1
else:
return compare_gamess_style(item1[:-1],item2[:-1])
def expend_and_order_sym(str_):
#Expend
for i,c in enumerate(str_):
try:
n = int(c)
except ValueError:
pass
else:
str_ = str_[:i-1] + str_[i-1]*n + str_[i+1:]
#Order by frequency
return "".join(sorted(str_,key=str_.count,reverse=True))
def get_nb_permutation(str_):
l = len(str_)-1
if l==0:
return 1
else:
return 2*(2*l + 1)
## We will order the symetry
l_sym_without_header = sym_raw.split("\n")[3:-2]
l_l_sym = [i.split() for i in l_sym_without_header]
for l in l_l_sym:
l[2] = expend_and_order_sym(l[2])
l_l_sym_iter = iter(l_l_sym)
for i,l in enumerate(l_l_sym_iter):
n = get_nb_permutation(l[2])
if n !=1:
l_l_sym[i:i+n] = sorted(l_l_sym[i:i+n],key=lambda x : x[2], cmp=compare_gamess_style)
for next_ in range(n-1):
next(l_l_sym_iter)
#Is orderd now
l_block = mo_raw.split("\n\n")[5:-1]
l_block_format=[]
print ""
print "BEGIN_MO"
for block in l_block:
print ""
l_ligne = block.split("\n")
print l_ligne.pop(0)
for l in l_l_sym:
i = int(l[0]) - 1
i_a = int(l[1]) - 1
sym = l[2]
print l_label[i_a],sym,l_ligne[i]
print "END_MO"
if do_pseudo:
print ""
print "BEGIN_PSEUDO"
klocmax = ezfio.get_pseudo_pseudo_klocmax()
kmax = ezfio.get_pseudo_pseudo_kmax()
lmax = ezfio.get_pseudo_pseudo_lmax()
n_k = ezfio.get_pseudo_pseudo_n_k()
v_k = ezfio.get_pseudo_pseudo_v_k()
dz_k = ezfio.get_pseudo_pseudo_dz_k()
n_kl = ezfio.get_pseudo_pseudo_n_kl()
v_kl = ezfio.get_pseudo_pseudo_v_kl()
dz_kl = ezfio.get_pseudo_pseudo_dz_kl()
def list_to_string(l):
return " ".join(map(str,l))
for i,a in enumerate(l_label):
l_str = []
l_dump = []
for k in range(klocmax):
if v_k[k][i]:
l_ = list_to_string([v_k[k][i], n_k[k][i]+2, dz_k[k][i]])
l_dump.append(l_)
l_str.append(l_dump)
for l in range(lmax+1):
l_dump = []
for k in range(kmax):
if v_kl[l][k][i]:
l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i]+2, dz_kl[l][k][i]])
l_dump.append(l_)
if l_dump:
l_str.append(l_dump)
str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE"
print str_.format(a,i+1,zcore[i],len(l_str))
for i, l in enumerate(l_str):
str_ = "FOR L= {0} COEFF N ZETA"
print str_.format(len(l_str)-i-1)
for ii, ll in enumerate(l):
print " ",ii+1, ll
str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS."
print str_.format(sum(zcore))
print "END_PSEUDO"
print ""
print "BEGIN_DET"
print ""
print "mo_num", mo_num
print "det_num", n_det
print ""
psi_det = ezfio.get_determinants_psi_det()
psi_coef = ezfio.get_determinants_psi_coef()[0]
for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef,psi_det):
print c
for det in l_det_bit_alpha:
bin_det_raw = "{0:b}".format(det)[::-1]
bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw))
print bin_det
for det in l_det_bit_beta:
bin_det_raw = "{0:b}".format(det)[::-1]
bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw))
print bin_det
print ""
print "END_DET"

359
plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py Executable file
View File

@ -0,0 +1,359 @@
#!/usr/bin/python
print "#QP -> QMCPACK"
# ___
# | ._ o _|_
# _|_ | | | |_
#
from ezfio import ezfio
import sys
ezfio_path = sys.argv[1]
ezfio.set_file(ezfio_path)
do_pseudo = ezfio.get_pseudo_do_pseudo()
if do_pseudo:
print "do_pseudo True"
zcore = ezfio.get_pseudo_nucl_charge_remove()
else:
print "do_pseudo False"
try:
n_det = ezfio.get_determinants_n_det()
except IOError:
n_det = 1
if n_det == 1:
print "multi_det False"
else:
print "multi_det True"
#
# |\/| o _ _
# | | | _> (_
#
def list_to_string(l):
return " ".join(map(str, l))
ao_num = ezfio.get_ao_basis_ao_num()
print "ao_num", ao_num
mo_num = ezfio.get_mo_basis_mo_tot_num()
print "mo_num", mo_num
alpha = ezfio.get_electrons_elec_alpha_num()
beta = ezfio.get_electrons_elec_beta_num()
print "elec_alpha_num", alpha
print "elec_beta_num", beta
print "elec_tot_num", alpha + beta
print "spin_multiplicity", 2 * (alpha - beta) + 1
l_label = ezfio.get_nuclei_nucl_label()
l_charge = ezfio.get_nuclei_nucl_charge()
l_coord = ezfio.get_nuclei_nucl_coord()
l_coord_str = [list_to_string(i) for i in zip(*l_coord)]
print "nucl_num", len(l_label)
# _
# / _ _ ._ _|
# \_ (_) (_) | (_|
#
print "Atomic coord in Bohr"
for i, t in enumerate(zip(l_label, l_charge, l_coord_str)):
try:
l = (t[0], t[1] + zcore[i], t[2])
except NameError:
l = t
print list_to_string(l)
#
# Call externet process to get the sysmetry
#
import subprocess
process = subprocess.Popen(
['qp_print_basis', ezfio_path],
stdout=subprocess.PIPE)
out, err = process.communicate()
basis_raw, sym_raw, _= out.split("\n\n\n")
# _ __
# |_) _. _ o _ (_ _ _|_
# |_) (_| _> | _> __) (/_ |_
#
basis_without_header = "\n".join(basis_raw.split("\n")[11:])
import re
l_basis_raw = re.split('\n\s*\n', basis_without_header)
a_already_print = []
l_basis_clean = []
for i, (a,b) in enumerate(zip(l_label,l_basis_raw)):
if a not in a_already_print:
l_basis_clean.append(b.replace('Atom {0}'.format(i + 1), a))
a_already_print.append(a)
else:
continue
print "BEGIN_BASIS_SET\n"
print "\n\n".join(l_basis_clean)
print "END_BASIS_SET"
# _
# |\/| / \ _
# | | \_/ _>
#
#
# Function
#
def same_character(item1):
return item1 == item1[0] * len(item1)
def compare_gamess_style(item1, item2):
if len(item1) < len(item2):
return -1
elif len(item1) > len(item2):
return 1
elif same_character(item1) and same_character(item2):
if item1 < item2:
return -1
else:
return 1
elif same_character(item1) and not same_character(item2):
return -1
elif not same_character(item1) and same_character(item2):
return 1
else:
return compare_gamess_style(item1[:-1], item2[:-1])
def expend_sym_str(str_):
#Expend x2 -> xx
# yx2 -> xxy
for i, c in enumerate(str_):
try:
n = int(c)
except ValueError:
pass
else:
str_ = str_[:i - 1] + str_[i - 1] * n + str_[i + 1:]
#Order by frequency
return "".join(sorted(str_, key=str_.count, reverse=True))
def expend_sym_l(l_l_sym):
for l in l_l_sym:
l[2] = expend_sym_str(l[2])
return l_l_sym
def get_nb_permutation(str_):
l = len(str_) - 1
if l == 0:
return 1
else:
return 2 * (2 * l + 1)
def order_l_l_sym(l_l_sym):
n = 1
for i in range(len(l_l_sym)):
if n != 1:
n += -1
continue
l = l_l_sym[i]
n = get_nb_permutation(l[2])
l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n],
key=lambda x: x[2],
cmp=compare_gamess_style)
return l_l_sym
#==========================
# We will order the symetry
#==========================
l_sym_without_header = sym_raw.split("\n")[3:-2]
l_l_sym_raw = [i.split() for i in l_sym_without_header]
l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw)
l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym)
#========
#MO COEF
#========
def order_phase(mo_coef):
#Order
mo_coef_phase = []
import math
for i in mo_coef:
if abs(max(i)) > abs(min(i)):
sign_max = math.copysign(1, max(i))
else:
sign_max = math.copysign(1, min(i))
if sign_max == -1:
ii = [-1 * l for l in i]
else:
ii = i
mo_coef_phase.append(ii)
return mo_coef_phase
def chunked(l, chunks_size):
l_block = []
for i in l:
chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)]
l_block.append(chunks)
return l_block
def print_mo_coef(mo_coef_block, l_l_sym):
print ""
print "BEGIN_MO"
print ""
len_block_curent = 0
nb_block = len(mo_coef_block[0])
for i_block in range(0, nb_block):
a = [i[i_block] for i in mo_coef_block]
r_ = range(len_block_curent, len_block_curent + len(a[0]))
print " ".join([str(i + 1) for i in r_])
len_block_curent += len(a[0])
for l in l_l_sym:
i = int(l[0]) - 1
i_a = int(l[1]) - 1
sym = l[2]
print l_label[i_a], sym, " ".join('{: 3.8f}'.format(i)
for i in a[i])
if i_block != nb_block - 1:
print ""
else:
print "END_MO"
mo_coef = ezfio.get_mo_basis_mo_coef()
mo_coef_transp = zip(*mo_coef)
mo_coef_block = chunked(mo_coef_transp, 4)
print_mo_coef(mo_coef_block, l_l_sym_ordered)
# _
# |_) _ _ _| _
# | _> (/_ |_| (_| (_)
#
if do_pseudo:
print ""
print "BEGIN_PSEUDO"
klocmax = ezfio.get_pseudo_pseudo_klocmax()
kmax = ezfio.get_pseudo_pseudo_kmax()
lmax = ezfio.get_pseudo_pseudo_lmax()
n_k = ezfio.get_pseudo_pseudo_n_k()
v_k = ezfio.get_pseudo_pseudo_v_k()
dz_k = ezfio.get_pseudo_pseudo_dz_k()
n_kl = ezfio.get_pseudo_pseudo_n_kl()
v_kl = ezfio.get_pseudo_pseudo_v_kl()
dz_kl = ezfio.get_pseudo_pseudo_dz_kl()
for i, a in enumerate(l_label):
l_str = []
#Local
l_dump = []
for k in range(klocmax):
if v_k[k][i]:
l_ = list_to_string([v_k[k][i], n_k[k][i] + 2, dz_k[k][i]])
l_dump.append(l_)
l_str.append(l_dump)
#Non local
for l in range(lmax + 1):
l_dump = []
for k in range(kmax):
if v_kl[l][k][i]:
l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i] + 2,
dz_kl[l][k][i]])
l_dump.append(l_)
if l_dump:
l_str.append(l_dump)
str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE"
print str_.format(a, i + 1, int(zcore[i]), int(len(l_str) - 1))
for i, l in enumerate(l_str):
str_ = "FOR L= {0} COEFF N ZETA"
print str_.format(int(len(l_str) - i - 1))
for ii, ll in enumerate(l):
print " ", ii + 1, ll
str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS."
print str_.format(sum(zcore))
print "END_PSEUDO"
# _
# | \ _ _|_
# |_/ (/_ |_
#
print ""
print "BEGIN_DET"
print ""
print "mo_num", mo_num
print "det_num", n_det
print ""
psi_det = ezfio.get_determinants_psi_det()
psi_coef = ezfio.get_determinants_psi_coef()[0]
for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef, psi_det):
print c
bin_det = ""
for i,int_det in enumerate(l_det_bit_alpha):
bin_det_raw = "{0:b}".format(int_det)[::-1]
if mo_num - 64*(i+1) > 0:
bin_det += bin_det_raw + "0" * (64*(i+1) - len(bin_det_raw))
else:
bin_det += bin_det_raw + "0" * (mo_num-64*i - len(bin_det_raw))
print bin_det
bin_det = ""
for i,int_det in enumerate(l_det_bit_beta):
bin_det_raw = "{0:b}".format(int_det)[::-1]
if mo_num - 64*(i+1) > 0:
bin_det += bin_det_raw + "0" * (64*(i+1) - len(bin_det_raw))
else:
bin_det += bin_det_raw + "0" * (mo_num-64*i - len(bin_det_raw))
print bin_det
print ""
print "END_DET"

2
plugins/qmcpack/save_for_qmcpack.irp.f
View File

@ -15,6 +15,6 @@ program qmcpack
enddo
call ezfio_set_ao_basis_ao_coef(ao_coef)
call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
call system('$QP_ROOT/src/qmcpack/qp_convert_qmcpack_from_ezfio.py '//trim(ezfio_filename))
call system('$QP_ROOT/src/qmcpack/qp_convert_qmcpack_to_ezfio.py '//trim(ezfio_filename))
end

140
scripts/entanglement.py Executable file
View File

@ -0,0 +1,140 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import matplotlib.pyplot as plt
from matplotlib import lines
from pylab import *
if __name__ == '__main__':
inputfile = sys.argv[1]
with open(str(inputfile),'r') as f:
vec_s1 = [float(x) for x in f.readline().split()]
mat_I = np.array([[float(x) for x in ln.split()] for ln in f])
#mat_I = np.loadtxt(str(inputfile),skiprows=1)
print vec_s1
print mat_I
#enable the following lines if you dont have a working X display
#import matplotlib
#matplotlib.use('Agg')
#MUTUALI INFORMATION PLOT:
# edit the following line to select the orbitals...
zselect=False
if zselect:
select =[]
select.extend(range(34,132))
select = [x-1 for x in select]
nselect = len(select)
vec_s1_s = np.zeros(nselect)
mat_I_s = np.zeros((nselect,nselect))
print mat_I_s, vec_s1_s
for i,x in enumerate(select):
vec_s1_s[i]=vec_s1[x]
for i,x in enumerate(select):
for j,y in enumerate(select):
mat_I_s[i,j] = mat_I[x,y]
print vec_s1_s
print mat_I_s
vec_s1 = vec_s1_s
mat_I = mat_I_s
#end selection
N= len(mat_I)
print N
theta = np.zeros(N)
r=np.zeros(N)
labels=np.zeros(N)
area=np.zeros(N)
for i in range(N):
theta[i]=-2*pi/N*i+pi/2
r[i]=1.0
if zselect:
labels[i]=select[i]+1
else:
labels[i]=i+1
area[i]=vec_s1[i]*500
if (len(sys.argv) > 2):
for j in range(2,len(sys.argv)):
if (sys.argv[j] == '-c'):
coordfile = sys.argv[j+1]
with open(str(coordfile),'r') as f:
for i in range(N):
line = f.readline().split()
theta[i] = float(line[0]) +pi/2
r[i] = float(line[1])
print 'theta=',theta
print 'r=',r
ax = plt.subplot(111,polar=True)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(b=False)
c = plt.scatter(theta,r,c="Red",s=area)
if (len(sys.argv) > 2):
for j in range(2,len(sys.argv)):
print 'sys.argv',j,sys.argv[j]
if (sys.argv[j] == '-t'):
plt.title(sys.argv[j+1])
break
else:
plt.title('N = '+str(N/2)+', Results file: '+sys.argv[1])
#this is dummy:
#c1 = plt.scatter(theta,(r+0.05),c="red",s=0)
for i in range(N):
# plt.annotate(int(labels[i]),xy=(theta[i],(r[i]+0.2)),size='xx-large',)
plt.text(theta[i],(r[i]+0.25),int(labels[i]),size='x-small',ha='center',va='center')
for j in range(i,N):
x =[theta[i],theta[j]]
y =[r[i],r[j]]
#y =[1,1]
if mat_I[i,j] >= 0.1:
line1 = lines.Line2D(x, y, linewidth=3, color='black',linestyle='-', alpha=0.8,label='0.1')
ax.add_line(line1)
elif mat_I[i,j] >=0.01:
line2 = lines.Line2D(x, y, linewidth=3, color='green',linestyle='-', alpha=0.8,label='0.01')
ax.add_line(line2)
elif mat_I[i,j] >=0.001:
line3 = lines.Line2D(x, y, linewidth=3, color='gray',linestyle='-', alpha=0.6,label='0.001')
ax.add_line(line3)
params = {'legend.fontsize': 16,
'legend.linewidth': 2}
#need to check if the three types of lines really exist.
if 'line1' in locals() and 'line2' in locals() and 'line3' in locals():
ax.legend([line1, line2, line3],[line1.get_label(),line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line1' in locals() and 'line2' in locals():
ax.legend([line1, line2],[line1.get_label(),line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line1' in locals() and 'line3' in locals():
ax.legend([line1, line3],[line1.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line2' in locals() and 'line3' in locals():
ax.legend([line2, line3],[line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line1' in locals():
ax.legend([line1],[line1.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line2' in locals():
ax.legend([line3],[line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
elif 'line3' in locals():
ax.legend([line2],[line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
#probably not the best way to code it.
#for saving, enable matplotlib.use('Agg')
#plt.savefig(str(sys.argv[1])+'.eps')
plt.show()

14
scripts/generate_h_apply.py
View File

@ -24,6 +24,8 @@ skip
init_main
filter_integrals
filter2h2p
filter_only_1h1p_single
filter_only_1h1p_double
filterhole
filterparticle
do_double_excitations
@ -150,6 +152,18 @@ class H_apply(object):
self["filterparticle"] = """
if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle
"""
def filter_only_1h1p(self):
self["filter_only_1h1p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(hole).eqv..False.) cycle
"""
self["filter_only_1h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(key).eqv..False.) cycle
"""
def unset_skip(self):
self["skip"] = """
"""

2
scripts/module/module_handler.py
View File

@ -40,7 +40,7 @@ def is_plugin(path_module_rel):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
return os.path.isfile(fpath) and os.access(fpath, os.X_OK) and not fpath.endswith(".py")
def get_dict_child(l_root_abs=None):

7
scripts/module/qp_module.py
View File

@ -175,7 +175,11 @@ def main(arguments):
d_child = d_local.copy()
d_child.update(d_plugin)
l_name = arguments["<name>"]
normalize_case = {}
for name in d_local.keys() + d_plugin.keys():
normalize_case [ name.lower() ] = name
l_name = [ normalize_case[name.lower()] for name in arguments["<name>"] ]
for name in l_name:
if name in d_local:
@ -206,6 +210,7 @@ def main(arguments):
raise
print "[ OK ]"
print ""
print "You can now compile as usual"
print "`cd {0} ; ninja` for exemple".format(QP_ROOT)
print " or --in developement mode-- you can cd in a directory and compile here"

27
src/AO_Basis/aos.irp.f
View File

@ -146,3 +146,30 @@ integer function ao_power_index(nx,ny,nz)
ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
end
BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
&BEGIN_PROVIDER [ integer, ao_l_max ]
&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ]
implicit none
BEGIN_DOC
! ao_l = l value of the AO: a+b+c in x^a y^b z^c
END_DOC
integer :: i
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
ao_l_char(i) = l_to_charater(ao_l(i))
enddo
ao_l_max = maxval(ao_l)
END_PROVIDER
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
BEGIN_DOC
! character corresponding to the "L" value of an AO orbital
END_DOC
implicit none
l_to_charater(0)='S'
l_to_charater(1)='P'
l_to_charater(2)='D'
l_to_charater(3)='F'
l_to_charater(4)='G'
END_PROVIDER

48
src/AO_Basis/aos_value.irp.f Normal file
View File

@ -0,0 +1,48 @@
double precision function ao_value(i,r)
implicit none
BEGIN_DOC
! return the value of the ith ao at point r
END_DOC
double precision, intent(in) :: r(3)
integer, intent(in) :: i
integer :: m,num_ao
double precision :: center_ao(3)
double precision :: beta
integer :: power_ao(3)
num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3)
double precision :: accu,dx,dy,dz,r2
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1)
dy = dy**power_ao(2)
dz = dz**power_ao(3)
accu = 0.d0
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
enddo
ao_value = accu * dx * dy * dz
end
subroutine give_all_aos_at_r(r,aos_array)
implicit none
BEGIN_dOC
! gives the values of aos at a given point r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
integer :: i
double precision :: ao_value
do i = 1, ao_num
aos_array(i) = ao_value(i,r)
enddo
end

44
src/Bitmask/bitmask_cas_routines.irp.f
View File

@ -445,3 +445,47 @@ integer function number_of_particles_verbose(key_in)
+ popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) )
end
logical function is_a_1h1p(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: number_of_particles, number_of_holes
is_a_1h1p = .False.
if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.1)then
is_a_1h1p = .True.
endif
end
logical function is_a_1h(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: number_of_particles, number_of_holes
is_a_1h = .False.
if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.0)then
is_a_1h = .True.
endif
end
logical function is_a_1p(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: number_of_particles, number_of_holes
is_a_1p = .False.
if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.1)then
is_a_1p = .True.
endif
end
logical function is_a_2p(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: number_of_particles, number_of_holes
is_a_2p = .False.
if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.2)then
is_a_2p = .True.
endif
end

89
src/Bitmask/bitmasks.irp.f
View File

@ -289,7 +289,12 @@ END_PROVIDER
&BEGIN_PROVIDER [ integer, n_virt_orb ]
implicit none
BEGIN_DOC
! Bitmasks for the inactive orbitals that are excited in post CAS method
! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited
! in post CAS methods
! n_inact_orb : Number of inactive orbitals
! virt_bitmask : Bitmaks of vritual orbitals which are supposed to be recieve electrons
! in post CAS methods
! n_virt_orb : Number of virtual orbitals
END_DOC
logical :: exists
integer :: j,i
@ -327,8 +332,14 @@ END_PROVIDER
BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)]
BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)]
&BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)]
BEGIN_DOC
! list_inact : List of the inactive orbitals which are supposed to be doubly excited
! in post CAS methods
! list_virt : List of vritual orbitals which are supposed to be recieve electrons
! in post CAS methods
END_DOC
implicit none
integer :: occ_inact(N_int*bit_kind_size)
integer :: itest,i
@ -348,6 +359,21 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
! Reunion of the inactive, active and virtual bitmasks
END_DOC
integer :: i,j
do i = 1, N_int
reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1))
reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2))
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
@ -376,7 +402,7 @@ END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
! Reunion of the inactive, active and virtual bitmasks
! Bitmask of the core orbitals that are never excited in post CAS method
END_DOC
integer :: i,j
do i = 1, N_int
@ -385,6 +411,35 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, list_core, (n_core_orb)]
BEGIN_DOC
! List of the core orbitals that are never excited in post CAS method
END_DOC
implicit none
integer :: occ_core(N_int*bit_kind_size)
integer :: itest,i
occ_core = 0
call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int)
ASSERT(itest==n_core_orb)
do i = 1, n_core_orb
list_core(i) = occ_core(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, n_core_orb ]
implicit none
BEGIN_DOC
! Number of core orbitals that are never excited in post CAS method
END_DOC
logical :: exists
integer :: j,i
integer :: i_hole,i_part,i_gen
n_core_orb = 0
do j = 1, N_int
n_core_orb += popcnt(core_bitmask(j,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, i_bitmask_gen ]
@ -407,3 +462,31 @@ END_PROVIDER
unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2))
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, n_act_orb]
BEGIN_DOC
! number of active orbitals
END_DOC
implicit none
integer :: i,j
n_act_orb = 0
do i = 1, N_int
n_act_orb += popcnt(cas_bitmask(i,1,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, list_act, (n_act_orb)]
BEGIN_DOC
! list of active orbitals
END_DOC
implicit none
integer :: occ_act(N_int*bit_kind_size)
integer :: itest,i
occ_act = 0
call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int)
ASSERT(itest==n_act_orb)
do i = 1, n_act_orb
list_act(i) = occ_act(i)
enddo
END_PROVIDER

6
src/Determinants/EZFIO.cfg
View File

@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2
interface: ezfio,provider,ocaml
default: False
[diagonalize_s2]
type: logical
doc: Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values.
interface: ezfio,provider,ocaml
default: True
[threshold_davidson]
type: Threshold
doc: Thresholds of Davidson's algorithm

5
src/Determinants/H_apply.template.f
View File

@ -165,6 +165,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
endif
logical :: check_double_excitation
logical :: is_a_1h1p
logical :: b_cycle
check_double_excitation = .True.
iproc = iproc_in
@ -298,6 +299,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l)
$filter2h2p
$filter_only_1h1p_double
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -348,6 +350,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibset(key(k,ispin),l)
$filter2h2p
$filter_only_1h1p_double
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -418,6 +421,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
integer(bit_kind) :: key_mask(N_int, 2)
logical :: check_double_excitation
logical :: is_a_1h1p
key_mask(:,:) = 0_bit_kind
@ -490,6 +494,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
$filterparticle
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
$filter2h2p
$filter_only_1h1p_single
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)

5
src/Determinants/connected_to_ref.irp.f
View File

@ -95,12 +95,13 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint)
exit
endif
enddo
i += 1
if (i > N_det) then
if (i >= N_det) then
return
endif
i += 1
!DIR$ FORCEINLINE
do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or. &

51
src/Determinants/density_matrix.irp.f
View File

@ -206,3 +206,54 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
state_average_weight = 1.d0/dble(N_states)
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_spin_density_ao, (ao_num_align,ao_num) ]
BEGIN_DOC
! one body spin density matrix on the AO basis : rho_AO(alpha) - rho_AO(beta)
END_DOC
implicit none
integer :: i,j,k,l
double precision :: dm_mo
one_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_spin_density_mo(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_spin_density_ao(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ]
BEGIN_DOC
! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta)
END_DOC
implicit none
integer :: i,j,k,l
double precision :: dm_mo
one_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_dm_mo_alpha(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
enddo
enddo
enddo
enddo
END_PROVIDER

191
src/Determinants/diagonalize_CI.irp.f
View File

@ -36,17 +36,36 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
implicit none
double precision :: ovrlp,u_dot_v
integer :: i_good_state
integer, allocatable :: index_good_state_array(:)
logical, allocatable :: good_state_array(:)
double precision, allocatable :: s2_values_tmp(:)
integer :: i_other_state
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
integer :: i_state
double precision :: s2,e_0
integer :: i,j,k
double precision, allocatable :: s2_eigvalues(:)
double precision, allocatable :: e_array(:)
integer, allocatable :: iorder(:)
do j=1,N_states_diag
! Guess values for the "N_states_diag" states of the CI_eigenvectors
do j=1,min(N_states_diag,N_det)
do i=1,N_det
CI_eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
do j=N_det+1,N_states_diag
do i=1,N_det
CI_eigenvectors(i,j) = 0.d0
enddo
enddo
if (diag_algorithm == "Davidson") then
@ -59,36 +78,77 @@ END_PROVIDER
else if (diag_algorithm == "Lapack") then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
CI_electronic_energy(:) = 0.d0
do i=1,N_det
CI_eigenvectors(i,1) = eigenvectors(i,1)
enddo
integer :: i_state
double precision :: s2
if (s2_eig) then
i_state = 0
allocate (s2_eigvalues(N_det))
allocate(index_good_state_array(N_det),good_state_array(N_det))
good_state_array = .False.
do j=1,N_det
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
print*,'s2 = ',s2
s2_eigvalues(j) = s2
! Select at least n_states states with S^2 values closed to "expected_s2"
if(dabs(s2-expected_s2).le.0.3d0)then
i_state += 1
do i=1,N_det
CI_eigenvectors(i,i_state) = eigenvectors(i,j)
enddo
CI_electronic_energy(i_state) = eigenvalues(j)
CI_eigenvectors_s2(i_state) = s2
i_state +=1
index_good_state_array(i_state) = j
good_state_array(j) = .True.
endif
if (i_state.ge.N_states_diag) then
exit
if(i_state.eq.N_states) then
exit
endif
enddo
if(i_state .ne.0)then
! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
enddo
CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
enddo
i_other_state = 0
do j = 1, N_det
if(good_state_array(j))cycle
i_other_state +=1
if(i_state+i_other_state.gt.n_states_diag)then
exit
endif
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
do i=1,N_det
CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
enddo
CI_electronic_energy(i_state+i_other_state) = eigenvalues(j)
CI_eigenvectors_s2(i_state+i_other_state) = s2
enddo
deallocate(index_good_state_array,good_state_array)
else
print*,''
print*,'!!!!!!!! WARNING !!!!!!!!!'
print*,' Within the ',N_det,'determinants selected'
print*,' and the ',N_states_diag,'states requested'
print*,' We did not find any state with S^2 values close to ',expected_s2
print*,' We will then set the first N_states eigenvectors of the H matrix'
print*,' as the CI_eigenvectors'
print*,' You should consider more states and maybe ask for diagonalize_s2 to be .True. or just enlarge the CI space'
print*,''
do j=1,min(N_states_diag,N_det)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy(j) = eigenvalues(j)
CI_eigenvectors_s2(j) = s2_eigvalues(j)
enddo
endif
deallocate(s2_eigvalues)
else
do j=1,N_states_diag
! Select the "N_states_diag" states of lowest energy
do j=1,min(N_det,N_states_diag)
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
@ -99,7 +159,100 @@ END_PROVIDER
endif
deallocate(eigenvectors,eigenvalues)
endif
if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
! Diagonalizing S^2 within the "n_states_diag" states found
allocate(s2_eigvalues(N_states_diag))
call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors,n_det,size(psi_det,3),size(CI_eigenvectors,1),min(n_states_diag,n_det),s2_eigvalues)
do j = 1, N_states_diag
do i = 1, N_det
psi_coef(i,j) = CI_eigenvectors(i,j)
enddo
enddo
if(s2_eig)then
! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value
! closer to the "expected_s2" set as input
allocate(index_good_state_array(N_det),good_state_array(N_det))
good_state_array = .False.
i_state = 0
do j = 1, N_states_diag
if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
good_state_array(j) = .True.
i_state +=1
index_good_state_array(i_state) = j
endif
enddo
! Sorting the i_state good states by energy
allocate(e_array(i_state),iorder(i_state))
do j = 1, i_state
do i = 1, N_det
CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(j))
enddo
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
CI_electronic_energy(j) = e_0
e_array(j) = e_0
iorder(j) = j
enddo
call dsort(e_array,iorder,i_state)
do j = 1, i_state
CI_electronic_energy(j) = e_array(j)
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(iorder(j)))
do i = 1, N_det
CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(iorder(j)))
enddo
! call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
! print*,'e = ',CI_electronic_energy(j)
! print*,'<e> = ',e_0
! call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),s2)
! print*,'s^2 = ',CI_eigenvectors_s2(j)
! print*,'<s^2>= ',s2
enddo
deallocate(e_array,iorder)
! Then setting the other states without any specific energy order
i_other_state = 0
do j = 1, N_states_diag
if(good_state_array(j))cycle
i_other_state +=1
do i = 1, N_det
CI_eigenvectors(i,i_state + i_other_state) = psi_coef(i,j)
enddo
CI_eigenvectors_s2(i_state + i_other_state) = s2_eigvalues(j)
call u0_H_u_0(e_0,CI_eigenvectors(1,i_state + i_other_state),n_det,psi_det,N_int)
CI_electronic_energy(i_state + i_other_state) = e_0
enddo
deallocate(index_good_state_array,good_state_array)
else
! Sorting the N_states_diag by energy, whatever the S^2 value is
allocate(e_array(n_states_diag),iorder(n_states_diag))
do j = 1, N_states_diag
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
e_array(j) = e_0
iorder(j) = j
enddo
call dsort(e_array,iorder,n_states_diag)
do j = 1, N_states_diag
CI_electronic_energy(j) = e_array(j)
do i = 1, N_det
CI_eigenvectors(i,j) = psi_coef(i,iorder(j))
enddo
CI_eigenvectors_s2(j) = s2_eigvalues(iorder(j))
enddo
deallocate(e_array,iorder)
endif
deallocate(s2_eigvalues)
endif
END_PROVIDER
subroutine diagonalize_CI

171
src/Determinants/s2.irp.f
View File

@ -214,4 +214,175 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
deallocate (shortcut, sort_idx, sorted, version)
end
subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys)
integer, intent(in) :: n,nmax_coefs,nmax_keys,nstates
double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates)
double precision, intent(out) :: s2(nstates,nstates)
double precision :: s2_tmp,accu
integer :: i,j,l,jj,ll,kk
integer, allocatable :: idx(:)
double precision, allocatable :: tmp(:,:)
BEGIN_DOC
! returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states
! psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j)
END_DOC
s2 = 0.d0
do ll = 1, nstates
do jj = 1, nstates
accu = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (i,j,kk,idx,tmp,s2_tmp) &
!$OMP SHARED (ll,jj,psi_keys_tmp,psi_coefs_tmp,N_int,n,nstates) &
!$OMP REDUCTION(+:accu)
allocate(idx(0:n))
!$OMP DO SCHEDULE(dynamic)
do i = 1, n
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(i,jj)
call filter_connected(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx)
do kk=1,idx(0)
j = idx(kk)
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(j,jj) + psi_coefs_tmp(i,jj) * s2_tmp * psi_coefs_tmp(j,ll)
enddo
enddo
!$OMP END DO NOWAIT
deallocate(idx)
!$OMP BARRIER
!$OMP END PARALLEL
s2(ll,jj) += accu
enddo
enddo
do i = 1, nstates
do j =i+1,nstates
accu = 0.5d0 * (s2(i,j) + s2(j,i))
s2(i,j) = accu
s2(j,i) = accu
enddo
enddo
end
subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nmax_coefs,nstates,s2_eigvalues)
BEGIN_DOC
! You enter with nstates vectors in psi_coefs_inout that may be coupled by S^2
! The subroutine diagonalize the S^2 operator in the basis of these states.
! The vectors that you obtain in output are no more coupled by S^2,
! which does not necessary mean that they are eigenfunction of S^2.
! n,nmax,nstates = number of determinants, physical dimension of the arrays and number of states
! keys_tmp = array of integer(bit_kind) that represents the determinants
! psi_coefs(i,j) = coeff of the ith determinant in the jth state
! VECTORS ARE SUPPOSED TO BE ORTHONORMAL IN INPUT
END_DOC
implicit none
use bitmasks
integer, intent(in) :: n,nmax_keys,nmax_coefs,nstates
integer(bit_kind), intent(in) :: keys_tmp(N_int,2,nmax_keys)
double precision, intent(inout) :: psi_coefs_inout(nmax_coefs,nstates)
!integer, intent(in) :: ndets_real,ndets_keys,ndets_coefs,nstates
!integer(bit_kind), intent(in) :: keys_tmp(N_int,2,ndets_keys)
!double precision, intent(inout) :: psi_coefs_inout(ndets_coefs,nstates)
double precision, intent(out) :: s2_eigvalues(nstates)
double precision,allocatable :: s2(:,:),overlap(:,:)
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
integer :: i,j,k
double precision, allocatable :: psi_coefs_tmp(:,:)
double precision :: accu,coef_contract
double precision :: u_dot_u,u_dot_v
print*,''
print*,'*********************************************************************'
print*,'Cleaning the various vectors by diagonalization of the S^2 matrix ...'
print*,''
print*,'nstates = ',nstates
allocate(s2(nstates,nstates),overlap(nstates,nstates))
do i = 1, nstates
overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
do j = i+1, nstates
overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
overlap(j,i) = overlap(i,j)
enddo
enddo
print*,'Overlap matrix in the basis of the states considered'
do i = 1, nstates
write(*,'(10(F16.10,X))')overlap(i,:)
enddo
call ortho_lowdin(overlap,size(overlap,1),nstates,psi_coefs_inout,size(psi_coefs_inout,1),n)
print*,'passed ortho'
do i = 1, nstates
overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
do j = i+1, nstates
overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
overlap(j,i) = overlap(i,j)
enddo
enddo
print*,'Overlap matrix in the basis of the Lowdin orthonormalized states '
do i = 1, nstates
write(*,'(10(F16.10,X))')overlap(i,:)
enddo
call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
print*,'S^2 matrix in the basis of the states considered'
double precision :: accu_precision_diag,accu_precision_of_diag
accu_precision_diag = 0.d0
accu_precision_of_diag = 0.d0
do i = 1, nstates
do j = i+1, nstates
if( ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then
s2(i,j) = 0.d0
s2(j,i) = 0.d0
endif
enddo
enddo
do i = 1, nstates
write(*,'(10(F10.6,X))')s2(i,:)
enddo
print*,'Diagonalizing the S^2 matrix'
allocate(eigvalues(nstates),eigvectors(nstates,nstates))
call lapack_diagd(eigvalues,eigvectors,s2,nstates,nstates)
print*,'Eigenvalues of s^2'
do i = 1, nstates
print*,'s2 = ',eigvalues(i)
s2_eigvalues(i) = eigvalues(i)
enddo
print*,'Building the eigenvectors of the S^2 matrix'
allocate(psi_coefs_tmp(nmax_coefs,nstates))
psi_coefs_tmp = 0.d0
do j = 1, nstates
do k = 1, nstates
coef_contract = eigvectors(k,j) ! <phi_k|Psi_j>
do i = 1, n_det
psi_coefs_tmp(i,j) += psi_coefs_inout(i,k) * coef_contract
enddo
enddo
enddo
do j = 1, nstates
accu = 0.d0
do i = 1, n_det
accu += psi_coefs_tmp(i,j) * psi_coefs_tmp(i,j)
enddo
print*,'Norm of vector = ',accu
accu = 1.d0/dsqrt(accu)
do i = 1, n_det
psi_coefs_inout(i,j) = psi_coefs_tmp(i,j) * accu
enddo
enddo
!call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
!print*,'S^2 matrix in the basis of the NEW states considered'
!do i = 1, nstates
! write(*,'(10(F16.10,X))')s2(i,:)
!enddo
deallocate(s2,eigvalues,eigvectors,psi_coefs_tmp,overlap)
end

27
src/Determinants/slater_rules.irp.f
View File

@ -1507,6 +1507,33 @@ subroutine get_occ_from_key(key,occ,Nint)
end
subroutine u0_H_u_0(e_0,u_0,n,keys_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint
double precision, intent(out) :: e_0
double precision, intent(in) :: u_0(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: H_jj(n)
double precision :: v_0(n)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
do i = 1, n
H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
enddo
call H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
e_0 = u_dot_v(v_0,u_0,n)/u_dot_u(u_0,n)
end
subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
use bitmasks
implicit none

283
src/Determinants/usefull_for_ovb.irp.f Normal file
View File

@ -0,0 +1,283 @@
integer function n_open_shell(det_in,nint)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: det_in(nint,2),nint
integer :: i
n_open_shell = 0
do i=1,Nint
n_open_shell += popcnt(iand(xor(det_in(i,1),det_in(i,2)),det_in(i,1)))
enddo
end
integer function n_closed_shell(det_in,nint)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: det_in(nint,2),nint
integer :: i
n_closed_shell = 0
do i=1,Nint
n_closed_shell += popcnt(iand(det_in(i,1),det_in(i,2)))
enddo
end
integer function n_closed_shell_cas(det_in,nint)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: det_in(nint,2),nint
integer(bit_kind) :: det_tmp(nint,2)
integer :: i
n_closed_shell_cas = 0
do i=1,Nint
det_tmp(i,1) = xor(det_in(i,1),reunion_of_core_inact_bitmask(i,1))
det_tmp(i,2) = xor(det_in(i,2),reunion_of_core_inact_bitmask(i,2))
enddo
!call debug_det(det_tmp,nint)
do i=1,Nint
n_closed_shell_cas += popcnt(iand(det_tmp(i,1),det_tmp(i,2)))
enddo
end
subroutine doubly_occ_empty_in_couple(det_in,n_couples,couples,couples_out)
implicit none
use bitmasks
integer, intent(in) :: n_couples,couples(n_couples,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: couples_out(0:n_couples)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
couples_out(0) = .False.
do i = 1, n_couples
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
! det_tmp is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 1)then
couples_out(i) = .False.
cycle
endif
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
couples_out(i) = .False.
else
couples_out(i) = .True.
couples_out(0) = .True.
endif
enddo
end
subroutine give_index_of_doubly_occ_in_active_space(det_in,doubly_occupied_array)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: det_in(N_int,2)
logical, intent(out) :: doubly_occupied_array(n_act_orb)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
do i = 1, n_act_orb
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
i_bite = list_act(i)
call set_bit_to_integer(i_bite,det_tmp_bis,N_int) ! act orb
! det_tmp is zero except for the orbital "ith" active orbital
integer :: i_count,i_bite
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
doubly_occupied_array(i) = .False.
else
doubly_occupied_array(i) = .True.
endif
enddo
end
subroutine doubly_occ_empty_in_couple_and_no_hund_elsewhere(det_in,n_couple_no_hund,couple_ion,couple_no_hund,is_ok)
implicit none
use bitmasks
integer, intent(in) :: n_couple_no_hund,couple_ion(2),couple_no_hund(n_couple_no_hund,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: is_ok
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
is_ok = .False.
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couple_ion(1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couple_ion(2),det_tmp_bis,N_int) ! second orb
! det_tmp is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 1)then
is_ok = .False.
return
endif
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
is_ok = .False.
return
else
do i = 1, n_couple_no_hund
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couple_no_hund (i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couple_no_hund (i,2),det_tmp_bis,N_int) ! second orb
! det_tmp_bis is zero except for the two orbitals of the couple
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 2)then
is_ok = .False.
return
endif
! test if orbital there are one alpha and one beta
integer :: i_count_alpha,i_count_beta
i_count_alpha = 0
i_count_beta = 0
do j = 1, N_int
i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
enddo
if(i_count_alpha==1.and.i_count_beta==1)then
is_ok = .True.
else
is_ok = .False.
return
endif
enddo
is_ok = .True.
endif
end
subroutine neutral_no_hund_in_couple(det_in,n_couples,couples,couples_out)
implicit none
use bitmasks
integer, intent(in) :: n_couples,couples(n_couples,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: couples_out(0:n_couples)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
couples_out(0) = .True.
do i = 1, n_couples
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
! det_tmp_bis is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 2)then
couples_out(i) = .False.
cycle
endif
! test if orbital there are one alpha and one beta
integer :: i_count_alpha,i_count_beta
i_count_alpha = 0
i_count_beta = 0
do j = 1, N_int
i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
enddo
if(i_count_alpha==1.and.i_count_beta==1)then
couples_out(i) = .True.
else
couples_out(i) = .False.
endif
enddo
do i = 1, n_couples
if(.not.couples_out(i))then
couples_out(0) = .False.
endif
enddo
end

3
src/Integrals_Bielec/.gitignore vendored
View File

@ -17,5 +17,4 @@ ZMQ
ezfio_interface.irp.f
irpf90.make
irpf90_entities
tags
test_integrals
tags

2
src/Integrals_Monoelec/.gitignore vendored
View File

@ -12,9 +12,7 @@ Makefile.depend
Nuclei
Pseudo
Utils
check_orthonormality
ezfio_interface.irp.f
irpf90.make
irpf90_entities
save_ortho_mos
tags

5
src/MOGuess/.gitignore vendored
View File

@ -4,7 +4,6 @@
AO_Basis
Electrons
Ezfio_files
H_CORE_guess
IRPF90_man
IRPF90_temp
Integrals_Monoelec
@ -15,8 +14,6 @@ Nuclei
Pseudo
Utils
ezfio_interface.irp.f
guess_overlap
irpf90.make
irpf90_entities
tags
truncate_mos
tags

150
src/MO_Basis/ao_ortho_canonical.irp.f
View File

@ -1,3 +1,87 @@
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ]
implicit none
BEGIN_DOC
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
END_DOC
integer :: i
integer, external :: ao_power_index
integer :: ibegin,j,k
ao_cart_to_sphe_coef(:,:) = 0.d0
! Assume order provided by ao_power_index
i = 1
ao_cart_to_sphe_num = 0
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_cart_to_sphe_num += 1
ao_cart_to_sphe_coef(i,ao_cart_to_sphe_num) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef(i+j-1,ao_cart_to_sphe_num+k) = cart_to_sphe_$SHELL(j,k)
enddo
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_num += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe'
end select
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ]
implicit none
BEGIN_DOC
! AO overlap matrix in the spherical basis set
END_DOC
double precision, allocatable :: S(:,:)
allocate (S(ao_cart_to_sphe_num,ao_num))
call dgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
ao_overlap,size(ao_overlap,1), 0.d0, &
S, size(S,1))
call dgemm('N','N',ao_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, 1.d0, &
S, size(S,1), &
ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, &
ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1))
deallocate(S)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_num) ]
implicit none
BEGIN_DOC
! AO_cart_to_sphe_coef^(-1)
END_DOC
call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, &
ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ]
implicit none
@ -8,47 +92,39 @@
END_DOC
integer :: i
ao_ortho_canonical_coef(:,:) = 0.d0
do i=1,ao_num
ao_ortho_canonical_coef(i,i) = 1.d0
enddo
if (ao_cartesian) then
do i=1,ao_num
ao_ortho_canonical_coef(i,i) = 1.d0
enddo
ao_ortho_canonical_num = ao_num
call ortho_canonical(ao_overlap,size(ao_overlap,1), &
ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), &
ao_ortho_canonical_num)
else
integer, external :: ao_power_index
integer :: ibegin,j,k
! Assume order provided by ao_power_index
i = 1
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_ortho_canonical_coef(i,i) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k)
enddo
enddo
i += size(cart_to_sphe_$SHELL,1)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in sphe_to_cart'
end select
double precision, allocatable :: S(:,:)
allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num))
S = 0.d0
do i=1,ao_cart_to_sphe_num
S(i,i) = 1.d0
enddo
ao_ortho_canonical_num = ao_cart_to_sphe_num
call ortho_canonical (ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num)
call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, &
ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), &
S, size(S,1), &
0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1))
deallocate(S)
endif
call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num)
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]

1
src/Utils/LinearAlgebra.irp.f
View File

@ -88,7 +88,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
endif
enddo
do i=m+1,n
print *, D(i)
D(i) = 0.d0
enddo

4
src/Utils/map_module.f90
View File

@ -69,7 +69,6 @@ subroutine cache_map_init(map,sze)
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_map_size_kind) :: sze
call omp_init_lock(map%lock)
call omp_set_lock(map%lock)
map%n_elements = 0_8
map%map_size = 0_8
@ -101,6 +100,9 @@ subroutine map_init(map,keymax)
stop 5
endif
sze = 2
do i=0_8,map%map_size
call omp_init_lock(map%map(i)%lock)
enddo
!$OMP PARALLEL DEFAULT(NONE) SHARED(map,sze) PRIVATE(i)
!$OMP DO SCHEDULE(STATIC,512)
do i=0_8,map%map_size

43
tests/bats/qp.bats
View File

@ -79,11 +79,39 @@ function run_FCI() {
eq $energy_pt2 $4 $thresh
}
function run_all_1h_1p() {
thresh=1.e-6
test_exe all_1h_1p || skip
ezfio set_file $1
ezfio set determinants n_det_max $2
ezfio set perturbation pt2_max $3
ezfio set determinants threshold_davidson 1.e-10
qp_run all_1h_1p $1 | tee $1.F1h1p.out
energy="$(ezfio get all_singles energy)"
eq $energy $4 $thresh
}
# ___
# | _ _ _|_
# | (/_ _> |_
#
#=== DHNO
@test "init DHNO chipman-dzp" {
run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio
}
@test "SCF DHNO chipman-dzp" {
run_HF dhno.ezfio -130.4278777822
}
@test "all_1h_1p DHNO chipman-dzp" {
qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio
run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202
}
#=== HBO
@test "init HBO STO-3G" {
run_init HBO.xyz "-b STO-3G" hbo.ezfio
@ -100,11 +128,11 @@ function run_FCI() {
}
@test "SCF H2O cc-pVDZ" {
run_HF h2o.ezfio -76.0270218692377
run_HF h2o.ezfio -0.760270218692179E+02
}
@test "FCI H2O cc-pVDZ" {
run_FCI h2o.ezfio 10000 -76.2382562245107 -76.2433933430971
run_FCI h2o.ezfio 10000 -0.762382562429778E+02 -0.762433933485226E+02
}
@test "CAS_SD H2O cc-pVDZ" {
@ -113,10 +141,10 @@ function run_FCI() {
ezfio set_file $INPUT
ezfio set perturbation do_pt2_end False
ezfio set determinants n_det_max 1000
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-25]"
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]"
qp_run cas_sd_selected $INPUT
energy="$(ezfio get cas_sd energy)"
eq $energy -76.2219816183745 1.E-6
eq $energy -0.762219854008117E+02 1.E-5
}
@test "MRCC H2O cc-pVDZ" {
@ -128,7 +156,8 @@ function run_FCI() {
ezfio set determinants read_wf True
qp_run mrcc_cassd $INPUT
energy="$(ezfio get mrcc_cassd energy)"
eq $energy -76.2313853324571 1.E-3
eq $energy -0.762303253805911E+02 1.E-3
}
@ -138,11 +167,11 @@ function run_FCI() {
}
@test "SCF H2O VDZ pseudo" {
run_HF h2o_pseudo.ezfio -16.9483703904632
run_HF h2o_pseudo.ezfio -0.169483703904991E+02
}
@test "FCI H2O VDZ pseudo" {
run_FCI h2o_pseudo.ezfio 2000 -17.1550015533975 -17.1645044211228
run_FCI h2o_pseudo.ezfio 2000 -0.171550015498807E+02 -0.171645044185009E+02
}
#=== Convert

7
tests/input/dhno.xyz Normal file
View File

@ -0,0 +1,7 @@
4
XYZ file: coordinates in Angstrom
H -0.877367 -1.047049 0.000000
N 0.000000 -0.544985 0.000000
O 0.000000 0.738624 0.000000
H 0.877367 -1.047049 0.000000

6
tests/run_tests.sh Executable file
View File

@ -0,0 +1,6 @@
#!/bin/bash
rm -rf work
exec bats bats/qp.bats