LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-31 08:36:02 +01:00
Code Issues Releases Wiki Activity

Changed Primitive.ml into GaussianPrimitive.ml

Browse Source
This commit is contained in:
Anthony Scemama 2017-05-31 23:47:18 +02:00
parent 840fe12c9d
commit 673747d530
17 changed files with 294 additions and 298 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
ocaml/Primitive.ml → ocaml/GaussianPrimitive.ml
View File

@ -1,5 +1,5 @@
open Qptypes;;
open Core.Std;;
open Qptypes
open Core.Std
type t =
{ sym : Symmetry.t ;
@ -11,8 +11,7 @@ let to_string p =
Printf.sprintf "(%s, %f)"
(Symmetry.to_string s)
(AO_expo.to_float e)
;;
let of_sym_expo s e =
{ sym=s ; expo=e}
;;

12
ocaml/Gto.ml
View File

@ -10,17 +10,17 @@ type fmt =
type t =
{ sym : Symmetry.t ;
lc : ((Primitive.t * AO_coef.t) list)
lc : ((GaussianPrimitive.t * AO_coef.t) list)
} with sexp
let of_prim_coef_list pc =
let (p,c) = List.hd_exn pc in
let sym = p.Primitive.sym in
let sym = p.GaussianPrimitive.sym in
let rec check = function
| [] -> `OK
| (p,c)::tl ->
if p.Primitive.sym <> sym then
if p.GaussianPrimitive.sym <> sym then
`Failed
else
check tl
@ -59,7 +59,7 @@ let read_one in_channel =
let coef = String.tr ~target:'D' ~replacement:'e' coef
in
let p =
Primitive.of_sym_expo sym
GaussianPrimitive.of_sym_expo sym
(AO_expo.of_float (Float.of_string expo) )
and c = AO_coef.of_float (Float.of_string coef) in
read_lines ( (p,c)::result) (i-1)
@ -80,7 +80,7 @@ let to_string_gamess { sym = sym ; lc = lc } =
let rec do_work accu i = function
| [] -> List.rev accu
| (p,c)::tail ->
let p = AO_expo.to_float p.Primitive.expo
let p = AO_expo.to_float p.GaussianPrimitive.expo
and c = AO_coef.to_float c
in
let result =
@ -100,7 +100,7 @@ let to_string_gaussian { sym = sym ; lc = lc } =
let rec do_work accu i = function
| [] -> List.rev accu
| (p,c)::tail ->
let p = AO_expo.to_float p.Primitive.expo
let p = AO_expo.to_float p.GaussianPrimitive.expo
and c = AO_coef.to_float c
in
let result =

6
ocaml/Gto.mli
View File

@ -6,12 +6,12 @@ type fmt =
type t =
{ sym : Symmetry.t ;
lc : (Primitive.t * Qptypes.AO_coef.t) list;
lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list;
} with sexp
(** Create from a list of Primitive.t * Qptypes.AO_coef.t *)
(** Create from a list of GaussianPrimitive.t * Qptypes.AO_coef.t *)
val of_prim_coef_list :
(Primitive.t * Qptypes.AO_coef.t) list -> t
(GaussianPrimitive.t * Qptypes.AO_coef.t) list -> t
(** Read from a file *)
val read_one : in_channel -> t

4
ocaml/Input_ao_basis.ml
View File

@ -112,8 +112,8 @@ end = struct
let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in
let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in
let prims = List.init ao_prim_num ~f:(fun j ->
let prim = { Primitive.sym = s ;
Primitive.expo = b.ao_expo.(ao_num*j+i)
let prim = { GaussianPrimitive.sym = s ;
GaussianPrimitive.expo = b.ao_expo.(ao_num*j+i)
}
in
let coef = b.ao_coef.(ao_num*j+i) in

22
ocaml/Pseudo.ml
View File

@ -2,7 +2,7 @@ open Core.Std
open Qptypes
module Primitive_local : sig
module GaussianPrimitive_local : sig
type t = {
expo : AO_expo.t ;
@ -29,7 +29,7 @@ end = struct
end
module Primitive_non_local : sig
module GaussianPrimitive_non_local : sig
type t = {
expo : AO_expo.t ;
@ -64,8 +64,8 @@ end
type t = {
element : Element.t ;
n_elec : Positive_int.t ;
local : (Primitive_local.t * AO_coef.t ) list ;
non_local : (Primitive_non_local.t * AO_coef.t ) list
local : (GaussianPrimitive_local.t * AO_coef.t ) list ;
non_local : (GaussianPrimitive_non_local.t * AO_coef.t ) list
} with sexp
let empty e =
@ -83,8 +83,8 @@ let to_string_local = function
( Printf.sprintf "%20s %8s %20s" "Coeff." "r^n" "Exp." ) ::
( List.map t ~f:(fun (l,c) -> Printf.sprintf "%20f %8d %20f"
(AO_coef.to_float c)
(R_power.to_int l.Primitive_local.r_power)
(AO_expo.to_float l.Primitive_local.expo)
(R_power.to_int l.GaussianPrimitive_local.r_power)
(AO_expo.to_float l.GaussianPrimitive_local.expo)
) )
|> String.concat ~sep:"\n"
@ -97,12 +97,12 @@ let to_string_non_local = function
( Printf.sprintf "%20s %8s %20s %8s" "Coeff." "r^n" "Exp." "Proj") ::
( List.map t ~f:(fun (l,c) ->
let p =
Positive_int.to_int l.Primitive_non_local.proj
Positive_int.to_int l.GaussianPrimitive_non_local.proj
in
Printf.sprintf "%20f %8d %20f |%d><%d|"
(AO_coef.to_float c)
(R_power.to_int l.Primitive_non_local.r_power)
(AO_expo.to_float l.Primitive_non_local.expo)
(R_power.to_int l.GaussianPrimitive_non_local.r_power)
(AO_expo.to_float l.GaussianPrimitive_non_local.expo)
p p
) )
|> String.concat ~sep:"\n"
@ -223,7 +223,7 @@ let read_element in_channel element =
let decode_local (pseudo,data) =
let decode_local_n n rest =
let result, rest =
loop Primitive_local.of_expo_r_power [] (Positive_int.to_int n,rest)
loop GaussianPrimitive_local.of_expo_r_power [] (Positive_int.to_int n,rest)
in
{ pseudo with local = result }, rest
in
@ -241,7 +241,7 @@ let read_element in_channel element =
let decode_non_local (pseudo,data) =
let decode_non_local_n proj n (pseudo,data) =
let result, rest =
loop (Primitive_non_local.of_proj_expo_r_power proj)
loop (GaussianPrimitive_non_local.of_proj_expo_r_power proj)
[] (Positive_int.to_int n, data)
in
{ pseudo with non_local = pseudo.non_local @ result }, rest

16
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -420,7 +420,7 @@ let run ?o b c d m p cart xyz_file =
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
Positive_int.to_int x.Pseudo.GaussianPrimitive_non_local.proj
in
if (x > accu) then x
else accu
@ -435,7 +435,7 @@ let run ?o b c d m p cart xyz_file =
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)
(Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i)
|> List.length )
|> List.fold ~init:0 ~f:(fun accu x ->
if accu > x then accu else x)
@ -458,8 +458,8 @@ let run ?o b c d m p cart xyz_file =
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
AO_expo.to_float y.Pseudo.GaussianPrimitive_local.expo,
R_power.to_int y.Pseudo.GaussianPrimitive_local.r_power
in
tmp_array_dz_k.(i).(j) <- y;
tmp_array_n_k.(i).(j) <- z;
@ -494,9 +494,9 @@ let run ?o b c d m p cart xyz_file =
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
Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj,
AO_expo.to_float y.Pseudo.GaussianPrimitive_non_local.expo,
R_power.to_int y.Pseudo.GaussianPrimitive_non_local.r_power
in
let i =
last_idx.(k)
@ -602,7 +602,7 @@ let run ?o b c d m p cart xyz_file =
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) )
prim.GaussianPrimitive.expo) )
end
in
let rec get_n n accu = function

23
ocaml/test_gto.ml
View File

@ -1,13 +1,13 @@
open Core.Std;;
open Qptypes;;
open Core.Std
open Qptypes
let test_prim () =
let p =
{ Primitive.sym = Symmetry.P ;
Primitive.expo = AO_expo.of_float 0.15} in
Primitive.to_string p
{ GaussianPrimitive.sym = Symmetry.P ;
GaussianPrimitive.expo = AO_expo.of_float 0.15} in
GaussianPrimitive.to_string p
|> print_string
;;
let test_gto_1 () =
let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
@ -26,24 +26,23 @@ let test_gto_1 () =
print_endline "gto3 = gto";
if (gto3 = gto3) then
print_endline "gto3 = gto3";
;;
let test_gto_2 () =
let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
ignore (input_line in_channel);
let basis = Basis.read in_channel (Nucl_number.of_int 1) in
List.iter basis ~f:(fun (x,n)-> Printf.printf "%d:%s\n" (Nucl_number.to_int n) (Gto.to_string x))
;;
let test_gto () =
let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
let basis = Basis.read_element in_channel (Nucl_number.of_int 1) Element.C in
List.iter basis ~f:(fun (x,n)-> Printf.printf "%d:%s\n" (Nucl_number.to_int n) (Gto.to_string x))
;;
let test_module () =
test_gto_1()
;;
test_module ();;
test_module ()

2
plugins/MRPT/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
MRPT_Utils Selectors_full Generators_full
MRPT_Utils Selectors_full Generators_full Psiref_CAS

158
plugins/MRPT_Utils/energies_cas.irp.f
View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
integer :: i
double precision :: energies(N_states)
do i = 1, N_states
call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i)
call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det_ref,psi_det_size,psi_det_size,N_states,i)
energy_cas_dyall(i) = energies(i)
print*, 'energy_cas_dyall(i)', energy_cas_dyall(i)
enddo
@ -15,7 +15,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
integer :: i
double precision :: energies(N_states)
do i = 1, N_states
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i)
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det_ref,psi_det_size,psi_det_size,N_states,i)
energy_cas_dyall_no_exchange(i) = energies(i)
print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i)
enddo
@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
double precision :: norm_out(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
use bitmasks
integer :: iorb
@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
orb = list_act(iorb)
hole_particle = 1
spin_exc = ispin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -53,8 +53,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
enddo
do state_target = 1,N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -72,7 +72,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb
integer :: state_target
@ -82,7 +82,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
orb = list_act(iorb)
hole_particle = -1
spin_exc = ispin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -93,8 +93,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -113,7 +113,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
orb_j = list_act(jorb)
hole_particle_j = 1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -139,10 +139,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
enddo
do state_target = 1 , N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -163,7 +163,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
@ -179,7 +179,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
orb_j = list_act(jorb)
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -189,10 +189,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
enddo
enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -213,7 +213,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
double precision :: energies(n_states)
@ -227,7 +227,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
orb_j = list_act(jorb)
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -238,14 +238,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
if(orb_i == orb_j .and. ispin .ne. jspin)then
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
else
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
endif
enddo
@ -268,7 +268,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -289,7 +289,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -301,12 +301,12 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
do state_target = 1, N_states
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -330,7 +330,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -351,7 +351,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -362,12 +362,12 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -391,7 +391,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -412,7 +412,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
orb_k = list_act(korb)
hole_particle_k = 1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -423,12 +423,12 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -452,7 +452,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -473,7 +473,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
@ -484,12 +484,12 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -515,7 +515,7 @@ END_PROVIDER
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb,i_ok
integer :: state_target
@ -541,10 +541,10 @@ END_PROVIDER
do state_target =1 , N_states
one_anhil_one_creat_inact_virt_norm(iorb,vorb,state_target,ispin) = 0.d0
enddo
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, N_int
psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,2,i) = psi_det(j,2,i)
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
if(i_ok.ne.1)then
@ -552,7 +552,7 @@ END_PROVIDER
call debug_det(psi_in_out,N_int)
print*, 'pb, i_ok ne 0 !!!'
endif
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
do j = 1, n_states
double precision :: coef,contrib
coef = psi_coef(i,j) !* psi_coef(i,j)
@ -585,7 +585,7 @@ END_PROVIDER
energies_alpha_beta(state_target, ispin) = - mo_bielec_integral_jj_exchange(orb_i,orb_v)
! energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
enddo
@ -620,7 +620,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: jorb,i_ok,aorb,orb_a
integer :: state_target
@ -645,10 +645,10 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
norm = 0.d0
norm_bis = 0.d0
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, N_int
psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,2,i) = psi_det(j,2,i)
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_a,ispin,i_ok)
if(i_ok.ne.1)then
@ -656,7 +656,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
psi_in_out_coef(i,j) = 0.d0
enddo
else
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
do j = 1, n_states
double precision :: coef,contrib
coef = psi_coef(i,j) !* psi_coef(i,j)
@ -688,7 +688,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
do state_target = 1, N_states
energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
enddo
@ -718,7 +718,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb,i_ok,aorb,orb_a
integer :: state_target
@ -743,10 +743,10 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
norm = 0.d0
norm_bis = 0.d0
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, N_int
psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,2,i) = psi_det(j,2,i)
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
call do_mono_excitation(psi_in_out(1,1,i),orb_a,orb_v,ispin,i_ok)
if(i_ok.ne.1)then
@ -754,7 +754,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
psi_in_out_coef(i,j) = 0.d0
enddo
else
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
do j = 1, n_states
double precision :: coef,contrib
coef = psi_coef(i,j) !* psi_coef(i,j)
@ -786,7 +786,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
do state_target = 1, N_states
energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
! print*, energies(state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
@ -826,7 +826,7 @@ END_PROVIDER
double precision, allocatable :: psi_in_out_coef(:,:)
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det+1,N_det+1))
allocate (eigenvectors(size(H_matrix,1),N_det+1))
allocate (eigenvalues(N_det+1))
@ -852,10 +852,10 @@ END_PROVIDER
- fock_virt_total_spin_trace(orb_v,j)
enddo
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, N_int
psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,2,i) = psi_det(j,2,i)
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
if(i_ok.ne.1)then
@ -865,7 +865,7 @@ END_PROVIDER
endif
interact_psi0(i) = 0.d0
do j = 1 , N_det
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,j),N_int,hij)
interact_psi0(i) += hij * psi_coef(j,1)
enddo
do j = 1, N_int
@ -975,7 +975,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:)
double precision, allocatable :: delta_e_det(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det+1,N_det+1))
allocate (eigenvectors(size(H_matrix,1),N_det+1))
allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det))
allocate (delta_e_det(N_det,N_det))
@ -1004,10 +1004,10 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
- fock_virt_total_spin_trace(orb_v,j)
enddo
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, N_int
psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,2,i) = psi_det(j,2,i)
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
if(i_ok.ne.1)then
@ -1017,8 +1017,8 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
endif
interact_psi0(i) = 0.d0
do j = 1 , N_det
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),delta_e_det(i,j))
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,j),N_int,hij)
call get_delta_e_dyall(psi_ref(1,1,j),psi_in_out(1,1,i),delta_e_det(i,j))
interact_cas(i,j) = hij
interact_psi0(i) += hij * psi_coef(j,1)
enddo

32
plugins/MRPT_Utils/mrpt_dress.irp.f
View File

@ -62,7 +62,7 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
call find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
if(N_tq > 0) then
call create_minilist(key_mask, psi_det, miniList, idx_miniList, N_det, N_minilist, Nint)
call create_minilist(key_mask, psi_ref, miniList, idx_miniList, N_det, N_minilist, Nint)
end if
@ -79,15 +79,15 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
phase_array =0.d0
do i = 1,idx_alpha(0)
index_i = idx_alpha(i)
call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha)
call i_h_j(tq(1,1,i_alpha),psi_ref(1,1,index_i),Nint,hialpha)
double precision :: coef_array(N_states)
do i_state = 1, N_states
coef_array(i_state) = psi_coef(index_i,i_state)
enddo
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
! call get_delta_e_dyall_general_mp(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e)
! call get_delta_e_dyall_general_mp(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e)
hij_array(index_i) = hialpha
call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
call get_excitation(psi_ref(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
! phase_array(index_i) = phase
do i_state = 1,N_states
delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
@ -100,12 +100,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
call omp_set_lock( psi_ref_bis_lock(index_i) )
do j = 1, idx_alpha(0)
index_j = idx_alpha(j)
! call get_excitation(psi_det(1,1,index_i),psi_det(1,1,index_i),exc,degree,phase,N_int)
! call get_excitation(psi_ref(1,1,index_i),psi_ref(1,1,index_i),exc,degree,phase,N_int)
! if(index_j.ne.index_i)then
! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then
! print*, phase_array(index_j) , phase_array(index_i) ,phase
! call debug_det(psi_det(1,1,index_i),N_int)
! call debug_det(psi_det(1,1,index_j),N_int)
! call debug_det(psi_ref(1,1,index_i),N_int)
! call debug_det(psi_ref(1,1,index_j),N_int)
! call debug_det(tq(1,1,i_alpha),N_int)
! stop
! endif
@ -123,14 +123,14 @@ end
BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
BEGIN_PROVIDER [ integer(bit_kind), gen_det_ref_sorted, (N_int,2,N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_ref_shortcut, (0:N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_ref_version, (N_int, N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_ref_idx, (N_det_generators,2) ]
gen_det_ref_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
gen_det_ref_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
call sort_dets_ab_v(gen_det_ref_sorted(:,:,:,1), gen_det_ref_idx(:,1), gen_det_ref_shortcut(0:,1), gen_det_ref_version(:,:,1), N_det_generators, N_int)
call sort_dets_ba_v(gen_det_ref_sorted(:,:,:,2), gen_det_ref_idx(:,2), gen_det_ref_shortcut(0:,2), gen_det_ref_version(:,:,2), N_det_generators, N_int)
END_PROVIDER

4
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -272,7 +272,7 @@ END_PROVIDER
allocate (s2_eigvalues(N_det))
allocate(index_good_state_array(N_det),good_state_array(N_det))
good_state_array = .False.
call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_ref,N_int,&
N_det,size(eigenvectors,1))
do j=1,N_det
! Select at least n_states states with S^2 values closed to "expected_s2"
@ -329,7 +329,7 @@ END_PROVIDER
deallocate(index_good_state_array,good_state_array)
deallocate(s2_eigvalues)
else
call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,&
call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_ref,N_int,&
min(N_det,N_states),size(eigenvectors,1))
! Select the "N_states" states of lowest energy
do j=1,min(N_det,N_states)

120
plugins/MRPT_Utils/new_way.irp.f
View File

@ -22,8 +22,8 @@ subroutine give_2h1p_contrib(matrix_2h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -45,7 +45,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
integer :: index_orb_act_mono(N_det,3)
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a)
@ -53,8 +53,8 @@ subroutine give_2h1p_contrib(matrix_2h1p)
do a = 1, n_act_orb ! First active
aorb = list_act(a)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -64,7 +64,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin
call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -81,7 +81,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_creat(a,jspin,istate) &
@ -109,7 +109,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a_{b} a^{\dagger}_a | Idet>
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_a
@ -150,7 +150,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
! you determine the interaction between the excited determinant and the other parent | Jdet >
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet >
! hja = < det_tmp | H | Jdet >
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,2) - active_int(borb,1) )
else
@ -216,8 +216,8 @@ subroutine give_1h2p_contrib(matrix_1h2p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -239,7 +239,7 @@ subroutine give_1h2p_contrib(matrix_1h2p)
integer :: index_orb_act_mono(N_det,3)
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb)
do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb)
@ -247,8 +247,8 @@ subroutine give_1h2p_contrib(matrix_1h2p)
aorb = list_act(a)
if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -258,7 +258,7 @@ subroutine give_1h2p_contrib(matrix_1h2p)
call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin
call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -280,7 +280,7 @@ subroutine give_1h2p_contrib(matrix_1h2p)
det_tmp(inint,2) = perturb_dets(inint,2,a,jspin,ispin)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_anhil(a,jspin,istate) &
@ -308,7 +308,7 @@ subroutine give_1h2p_contrib(matrix_1h2p)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a^{\dagger}_b a_{a} | Idet>
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a
@ -350,7 +350,7 @@ subroutine give_1h2p_contrib(matrix_1h2p)
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{borb,kspin} a_{iorb,ispin} | Jdet >
! hja = < det_tmp | H | Jdet >
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,1) - active_int(borb,2) )
else
@ -418,8 +418,8 @@ subroutine give_1h1p_contrib(matrix_1h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -430,20 +430,20 @@ subroutine give_1h1p_contrib(matrix_1h1p)
- fock_virt_total_spin_trace(rorb,j)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
do jdet = 1, idx(0)
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
double precision :: himono,delta_e(N_states),coef_mono(N_states)
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,himono)
do state_target = 1, N_states
! delta_e(state_target) = one_anhil_one_creat_inact_virt(i,r,state_target) + delta_e_inact_virt(state_target)
@ -451,7 +451,7 @@ subroutine give_1h1p_contrib(matrix_1h1p)
coef_mono(state_target) = himono / delta_e(state_target)
enddo
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
aorb = (exc(1,2,1)) !!! a^{\dagger}_a
@ -464,15 +464,15 @@ subroutine give_1h1p_contrib(matrix_1h1p)
jspin = 2
endif
call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
call get_excitation_degree(psi_ref(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
if(degree_scalar .ne. 2)then
print*, 'pb !!!'
print*, degree_scalar
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
call debug_det(det_tmp,N_int)
stop
endif
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(ispin == jspin )then
hij = -get_mo_bielec_integral(iorb,aorb,rorb,borb,mo_integrals_map) &
+ get_mo_bielec_integral(iorb,aorb,borb,rorb,mo_integrals_map)
@ -482,17 +482,17 @@ subroutine give_1h1p_contrib(matrix_1h1p)
hij = hij * phase
double precision :: hij_test
integer :: state_target
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test)
if(dabs(hij - hij_test).gt.1.d-10)then
print*, 'ahah pb !!'
print*, 'hij .ne. hij_test'
print*, hij,hij_test
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
call debug_det(det_tmp,N_int)
print*, ispin, jspin
print*,iorb,borb,rorb,aorb
print*, phase
call i_H_j_verbose(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
call i_H_j_verbose(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test)
stop
endif
do state_target = 1, N_states
@ -542,13 +542,13 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
double precision :: himono,delta_e(N_states),coef_mono(N_states)
integer :: state_target
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
do r = 1, n_virt_orb ! First virtual
@ -563,13 +563,13 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
- fock_virt_total_spin_trace(rorb,j)
enddo
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,himono)
do inint = 1, N_int
det_pert(inint,1,i,r,ispin) = det_tmp(inint,1)
det_pert(inint,2,i,r,ispin) = det_tmp(inint,2)
@ -630,7 +630,7 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
do jdet = 1, idx(0)
!
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
aorb = (exc(1,2,1)) !!! a^{\dagger}_a
@ -642,24 +642,24 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
jspin = 2
endif
call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
call get_excitation_degree(psi_ref(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
if(degree_scalar .ne. 2)then
print*, 'pb !!!'
print*, degree_scalar
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
call debug_det(det_tmp,N_int)
stop
endif
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
double precision :: hij_test
hij_test = 0.d0
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test)
do state_target = 1, N_states
matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
enddo
else
hij_test = 0.d0
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij_test)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hij_test)
do state_target = 1, N_states
matrix_1h1p(idet,idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
enddo
@ -701,13 +701,13 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
double precision :: himono,delta_e(N_states),coef_mono(N_states)
integer :: state_target
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
do i = 1, n_act_orb ! First active
iorb = list_act(i)
do r = 1, n_virt_orb ! First virtual
@ -721,8 +721,8 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p)
delta_e_act_virt(j) = - fock_virt_total_spin_trace(rorb,j)
enddo
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation active -- > virtual
call do_mono_excitation(det_tmp,iorb,rorb,ispin,i_ok)
@ -739,7 +739,7 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p)
enddo
cycle
endif
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,himono)
do inint = 1, N_int
det_pert(inint,1,i,r,ispin) = det_tmp(inint,1)
det_pert(inint,2,i,r,ispin) = det_tmp(inint,2)
@ -803,8 +803,8 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p)
enddo
do jdet = 1,N_det
double precision :: coef_array(N_states),hij_test
call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono)
call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,delta_e)
call i_H_j(det_tmp,psi_ref(1,1,jdet),N_int,himono)
call get_delta_e_dyall(psi_ref(1,1,jdet),det_tmp,delta_e)
do state_target = 1, N_states
! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1)
matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target)
@ -850,8 +850,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -862,7 +862,7 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
- fock_virt_total_spin_trace(rorb,j)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_double_alpha_beta(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_double_alpha_beta(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
do ispin = 1, 2
@ -872,8 +872,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
do b = 1, n_act_orb
borb = list_act(b)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation (i-->a)(ispin) + (b-->r)(other_spin(ispin))
integer :: i_ok,corb,dorb
@ -904,7 +904,7 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
pert_det(inint,2,a,b,ispin) = det_tmp(inint,2)
enddo
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hidouble)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hidouble)
do state_target = 1, N_states
delta_e(state_target) = one_anhil_one_creat(a,b,ispin,jspin,state_target) + delta_e_inact_virt(state_target)
pert_det_coef(a,b,ispin,state_target) = hidouble / delta_e(state_target)
@ -915,7 +915,7 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
enddo
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
integer :: c,d,state_target
integer(bit_kind) :: det_tmp_bis(N_int,2)
! excitation from I --> J
@ -935,8 +935,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
det_tmp_bis(inint,2) = pert_det(inint,2,c,d,2)
enddo
double precision :: hjdouble_1,hjdouble_2
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1)
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2)
call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1)
call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2)
do state_target = 1, N_states
matrix_1h1p(idx(jdet),idet,state_target) += (pert_det_coef(c,d,1,state_target) * hjdouble_1 + pert_det_coef(c,d,2,state_target) * hjdouble_2 )
enddo

80
plugins/MRPT_Utils/new_way_second_order_coef.irp.f
View File

@ -24,8 +24,8 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -50,9 +50,9 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
integer :: istate
do idet = 1, N_det
call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono_or_exchange(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
! if(idet == 81)then
! call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx)
! call get_excitation_degree_vector_mono_or_exchange_verbose(psi_ref(1,1,1),psi_ref(1,1,idet),degree,N_int,N_det,idx)
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
@ -61,8 +61,8 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
do a = 1, n_act_orb ! First active
aorb = list_act(a)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -72,7 +72,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin
call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -90,7 +90,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_creat(a,jspin,istate) &
@ -124,7 +124,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
@ -149,7 +149,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
else if(degree(jdet)==2)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b}
@ -196,7 +196,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! ! < idet | H | det_tmp > = phase * (ir|cv)
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -215,7 +215,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
stop
endif
! < jdet | H | det_tmp_bis > = phase * (ir|cv)
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -243,7 +243,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -260,7 +260,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
stop
endif
! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -296,7 +296,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
! | det_tmp > = a^{\dagger}_{aorb,beta} | Idet >
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
call get_double_excitation(det_tmp,psi_ref(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(aorb,1) - active_int(aorb,2))
else
@ -312,15 +312,15 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
else if(index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),4))then !! closed shell double excitation
else
call get_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,degree_scalar,phase,N_int)
call get_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,degree_scalar,phase,N_int)
integer :: h1,h2,p1,p2,s1,s2 , degree_scalar
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*, h1,p1,h2,p2,s1,s2
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idet),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
print*, idet,idx(jdet)
print*, 'pb !!!!!!!!!!!!!'
call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono_or_exchange_verbose(psi_ref(1,1,1),psi_ref(1,1,idet),degree,N_int,N_det,idx)
stop
endif
endif
@ -398,8 +398,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -430,7 +430,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
+ fock_core_inactive_total_spin_trace(iorb,istate)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono_or_exchange(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb)
do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb)
@ -443,8 +443,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
enddo
if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -454,7 +454,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin
call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -477,7 +477,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp(inint,2) = perturb_dets(inint,2,a,jspin,ispin)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
@ -501,7 +501,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
@ -526,7 +526,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
else if(degree(jdet)==2)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b}
@ -575,7 +575,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * (ir|cv)
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -590,7 +590,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = (fock_operator_local(aorb,borb,kspin) ) * phase
! < jdet | H | det_tmp_bis > = phase * (ir|cv)
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -617,7 +617,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -630,7 +630,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
! ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = fock_operator_local(aorb,borb,kspin) * phase
! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
@ -665,7 +665,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
call get_double_excitation(det_tmp,psi_ref(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(borb,1) - active_int(borb,2))
else
@ -674,8 +674,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then
call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),1,i_ok)
if(i_ok .ne. 1)then
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idet),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
print*, aorb, borb
call debug_det(det_tmp,N_int)
stop
@ -692,7 +692,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
stop
endif
hab = fock_operator_local(aorb,borb,1) * phase
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(borb,1) - active_int(borb,2))
else
@ -716,7 +716,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
call get_double_excitation(det_tmp,psi_ref(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(borb,1) - active_int(borb,2))
else
@ -725,8 +725,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then
call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),2,i_ok)
if(i_ok .ne. 1)then
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call debug_det(psi_ref(1,1,idet),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
print*, aorb, borb
call debug_det(det_tmp,N_int)
stop
@ -739,7 +739,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = fock_operator_local(aorb,borb,2) * phase
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(det_tmp_bis,psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(borb,1) - active_int(borb,2))
else

4
plugins/MRPT_Utils/psi_active_prov.irp.f
View File

@ -11,8 +11,8 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
!print*, 'psi_active '
do i = 1, N_det
do j = 1, N_int
psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1))
psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1))
psi_active(j,1,i) = iand(psi_ref(j,1,i),cas_bitmask(j,1,1))
psi_active(j,2,i) = iand(psi_ref(j,2,i),cas_bitmask(j,1,1))
enddo
enddo
END_PROVIDER

58
plugins/MRPT_Utils/second_order_new.irp.f
View File

@ -51,8 +51,8 @@ subroutine give_1h2p_new(matrix_1h2p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -79,7 +79,7 @@ subroutine give_1h2p_new(matrix_1h2p)
+ fock_core_inactive_total_spin_trace(iorb,istate)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono_or_exchange(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb)
do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb)
@ -90,8 +90,8 @@ subroutine give_1h2p_new(matrix_1h2p)
enddo
if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -101,7 +101,7 @@ subroutine give_1h2p_new(matrix_1h2p)
call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin
call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -116,7 +116,7 @@ subroutine give_1h2p_new(matrix_1h2p)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
@ -138,7 +138,7 @@ subroutine give_1h2p_new(matrix_1h2p)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | K_{ab} | Idet>
do jdet = 1, idx(0)
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
@ -163,7 +163,7 @@ subroutine give_1h2p_new(matrix_1h2p)
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
else if(degree(jdet)==2)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a ALPHA
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b} ALPHA
@ -209,13 +209,13 @@ subroutine give_1h2p_new(matrix_1h2p)
det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hia = phase * (active_int(aorb,1) - active_int(aorb,2) )
else
hia = phase * active_int(aorb,1)
endif
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,1) - active_int(borb,2) )
else
@ -254,7 +254,7 @@ subroutine give_1h2p_new(matrix_1h2p)
if(dabs(hia).le.1.d-12)cycle
if(dabs(hab).le.1.d-12)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(jspin == ispin)then
hjb = phase * (active_int(corb,1) - active_int(corb,2) )
else
@ -307,7 +307,7 @@ subroutine give_1h2p_new(matrix_1h2p)
hab = fock_operator_local(aorb,borb,1) * phase
if(dabs(hab).le.1.d-12)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(ispin == 2)then
hjb = phase * (active_int(aorb,1) - active_int(aorb,2) )
else if (ispin == 1)then
@ -341,7 +341,7 @@ subroutine give_1h2p_new(matrix_1h2p)
hab = fock_operator_local(aorb,borb,2) * phase
if(dabs(hab).le.1.d-12)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(ispin == 1)then
hjb = phase * (active_int(borb,1) - active_int(borb,2) )
else if (ispin == 2)then
@ -380,7 +380,7 @@ subroutine give_1h2p_new(matrix_1h2p)
hab = fock_operator_local(aorb,borb,1) * phase
if(dabs(hab).le.1.d-12)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(ispin == 2)then
hjb = phase * (active_int(borb,1) - active_int(borb,2) )
else if (ispin == 1)then
@ -415,7 +415,7 @@ subroutine give_1h2p_new(matrix_1h2p)
hab = fock_operator_local(aorb,borb,2) * phase
if(dabs(hab).le.1.d-12)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(ispin == 1)then
hjb = phase * (active_int(borb,1) - active_int(borb,2) )
else if (ispin == 2)then
@ -433,9 +433,9 @@ subroutine give_1h2p_new(matrix_1h2p)
else
! one should not fall in this case ...
call debug_det(psi_det(1,1,i),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call debug_det(psi_ref(1,1,i),N_int)
call debug_det(psi_ref(1,1,idx(jdet)),N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
call decode_exc(exc,2,h1,p1,h2,p2,s1,s2)
integer :: h1, p1, h2, p2, s1, s2
print*, h1, p1, h2, p2, s1, s2
@ -519,8 +519,8 @@ subroutine give_2h1p_new(matrix_2h1p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
@ -547,7 +547,7 @@ subroutine give_2h1p_new(matrix_2h1p)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a)
@ -555,8 +555,8 @@ subroutine give_2h1p_new(matrix_2h1p)
do a = 1, n_act_orb ! First active
aorb = list_act(a)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
@ -566,7 +566,7 @@ subroutine give_2h1p_new(matrix_2h1p)
call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin
call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
@ -580,7 +580,7 @@ subroutine give_2h1p_new(matrix_2h1p)
perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_creat(a,jspin,istate) + delta_e_inactive_virt(istate)
@ -602,7 +602,7 @@ subroutine give_2h1p_new(matrix_2h1p)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a_{b} a^{\dagger}_a | Idet>
do jdet = 1, idx(0)
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_part = list_act_reverse(exc(1,2,1)) ! a^{\dagger}_{aorb}
@ -658,7 +658,7 @@ subroutine give_2h1p_new(matrix_2h1p)
! you determine the interaction between the excited determinant and the other parent | Jdet >
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet >
! hja = < det_tmp | H | Jdet >
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,1) - active_int(borb,2) )
else
@ -698,7 +698,7 @@ subroutine give_2h1p_new(matrix_2h1p)
hab = fock_operator_local(borb,aorb,kspin) * phase
if(dabs(hab).le.1.d-10)cycle
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(jspin == ispin)then
hjb = phase * (active_int(corb,1) - active_int(corb,2) )
else

28
plugins/MRPT_Utils/second_order_new_2p.irp.f
View File

@ -50,8 +50,8 @@ subroutine give_2p_new(matrix_2p)
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
enddo
do v = 1, n_virt_orb ! First virtual
vorb = list_virt(v)
@ -82,8 +82,8 @@ subroutine give_2p_new(matrix_2p)
- fock_virt_total_spin_trace(vorb,istate)
enddo
do idet = 1, N_det
! call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
! call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (aorb,rorb)
do jspin = 1, 2 ! spin of the couple a-a^dagger (borb,vorb)
@ -108,8 +108,8 @@ subroutine give_2p_new(matrix_2p)
cycle ! condition not to double count
endif
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
det_tmp(inint,1) = psi_ref(inint,1,idet)
det_tmp(inint,2) = psi_ref(inint,2,idet)
enddo
! Do the excitation (aorb,ispin) --> (rorb,ispin)
call clear_bit_to_integer(aorb,det_tmp(1,ispin),N_int) ! hole in "aorb" of spin Ispin
@ -119,7 +119,7 @@ subroutine give_2p_new(matrix_2p)
call clear_bit_to_integer(borb,det_tmp(1,jspin),N_int) ! hole in "borb" of spin Jspin
call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
! Check if the excitation is possible or not on psi_ref(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,1)) + popcnt(det_tmp(inint,2))
@ -134,7 +134,7 @@ subroutine give_2p_new(matrix_2p)
perturb_dets(inint,2,a,b,ispin,jspin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,b,ispin,jspin) = phase
do istate = 1, N_states
@ -146,16 +146,16 @@ subroutine give_2p_new(matrix_2p)
else
perturb_dets_hij(a,b,ispin,jspin) = phase * active_int(a,b,1)
endif
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij)
call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hij)
if(hij.ne.perturb_dets_hij(a,b,ispin,jspin))then
print*, active_int(a,b,1) , active_int(b,a,1)
double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,idet),det_tmp,N_int,hij,hmono,hdouble)
call i_H_j_verbose(psi_ref(1,1,idet),det_tmp,N_int,hij,hmono,hdouble)
print*, 'pb !! hij.ne.perturb_dets_hij(a,b,ispin,jspin)'
print*, ispin,jspin
print*, aorb,rorb,borb,vorb
print*, hij,perturb_dets_hij(a,b,ispin,jspin)
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_ref(1,1,idet),N_int)
call debug_det(det_tmp,N_int)
stop
endif
@ -170,7 +170,7 @@ subroutine give_2p_new(matrix_2p)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | K_{ab} | Idet>
do jdet = 1, idx(0)
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
@ -195,7 +195,7 @@ subroutine give_2p_new(matrix_2p)
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
else if(degree(jdet)==2)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a ALPHA
@ -262,7 +262,7 @@ subroutine give_2p_new(matrix_2p)
do jdet = 1, idx(0)
! if(idx(jdet).gt.idet)cycle
do istate = 1, N_states
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij)
call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij)
matrix_2p(idx(jdet),idet,istate) += hij * perturb_dets_hij(a,b,ispin,jspin) * delta_e_inv(a,b,ispin,jspin,istate)
enddo
enddo ! jdet

16
plugins/mrcepa0/dressing.irp.f
View File

@ -99,6 +99,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
logical, external :: detEq, is_generable
!double precision, external :: get_dij, get_dij_index
if (perturbative_triples) then
PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
endif
leng = max(N_det_generators, N_det_non_ref)
@ -191,14 +194,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
end do
end if
if (perturbative_triples) then
double precision :: Delta_E_inv(N_states)
double precision, external :: diag_H_mat_elem
do i_state=1,N_states
Delta_E_inv(i_state) = 1.d0 / (psi_ref_energy_diagonalized(i_state) - diag_H_mat_elem(tq(1,1,i_alpha),N_int) )
enddo
endif
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
@ -263,9 +258,12 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
else if (perturbative_triples) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
hka = hij_cache(idx_alpha(k_sd))
do i_state=1,N_states
dka(i_state) = hka * Delta_E_inv(i_state)
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif