mirror of
https://github.com/LCPQ/quantum_package
synced 2024-06-25 22:52:15 +02:00
merged with toto
This commit is contained in:
Emmanuel Giner
commit
03c50ebd2f
|
|
@ -10,7 +10,7 @@
|
|||
#
|
||||
#
|
||||
[COMMON]
|
||||
FC : gfortran -ffree-line-length-none -I . -mavx -g
|
||||
FC : gfortran -ffree-line-length-none -I . -mavx -g
|
||||
LAPACK_LIB : -llapack -lblas
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32
|
||||
|
|
@ -35,7 +35,7 @@ OPENMP : 1 ; Append OpenMP flags
|
|||
# -ffast-math and the Fortran-specific
|
||||
# -fno-protect-parens and -fstack-arrays.
|
||||
[OPT]
|
||||
FCFLAGS : -Ofast
|
||||
FCFLAGS : -Ofast -march=native
|
||||
|
||||
# Profiling flags
|
||||
#################
|
||||
|
|
|
|||
|
|
@ -51,7 +51,7 @@ FCFLAGS : -Ofast
|
|||
# -g : Extra debugging information
|
||||
#
|
||||
[DEBUG]
|
||||
FCFLAGS : -g -msse4.2
|
||||
FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
|
||||
|
||||
# OpenMP flags
|
||||
#################
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@ module Determinants_by_hand : sig
|
|||
{ n_int : N_int_number.t;
|
||||
bit_kind : Bit_kind.t;
|
||||
n_det : Det_number.t;
|
||||
n_states : States_number.t;
|
||||
expected_s2 : Positive_float.t;
|
||||
psi_coef : Det_coef.t array;
|
||||
psi_det : Determinant.t array;
|
||||
|
|
@ -18,11 +19,14 @@ module Determinants_by_hand : sig
|
|||
val to_rst : t -> Rst_string.t
|
||||
val of_rst : Rst_string.t -> t option
|
||||
val read_n_int : unit -> N_int_number.t
|
||||
val update_ndet : Det_number.t -> unit
|
||||
val extract_state : States_number.t -> unit
|
||||
end = struct
|
||||
type t =
|
||||
{ n_int : N_int_number.t;
|
||||
bit_kind : Bit_kind.t;
|
||||
n_det : Det_number.t;
|
||||
n_states : States_number.t;
|
||||
expected_s2 : Positive_float.t;
|
||||
psi_coef : Det_coef.t array;
|
||||
psi_det : Determinant.t array;
|
||||
|
|
@ -129,12 +133,12 @@ end = struct
|
|||
|> Array.map ~f:Det_coef.of_float
|
||||
;;
|
||||
|
||||
let write_psi_coef ~n_det c =
|
||||
let write_psi_coef ~n_det ~n_states c =
|
||||
let n_det = Det_number.to_int n_det
|
||||
and c = Array.to_list c
|
||||
|> List.map ~f:Det_coef.to_float
|
||||
and n_states =
|
||||
read_n_states () |> States_number.to_int
|
||||
States_number.to_int n_states
|
||||
in
|
||||
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
|
||||
|> Ezfio.set_determinants_psi_coef
|
||||
|
|
@ -200,6 +204,7 @@ end = struct
|
|||
expected_s2 = read_expected_s2 () ;
|
||||
psi_coef = read_psi_coef () ;
|
||||
psi_det = read_psi_det () ;
|
||||
n_states = read_n_states () ;
|
||||
}
|
||||
else
|
||||
failwith "No molecular orbitals, so no determinants"
|
||||
|
|
@ -222,12 +227,14 @@ end = struct
|
|||
expected_s2 ;
|
||||
psi_coef ;
|
||||
psi_det ;
|
||||
n_states ;
|
||||
} =
|
||||
write_n_int n_int ;
|
||||
write_bit_kind bit_kind;
|
||||
write_n_det n_det;
|
||||
write_n_states n_states;
|
||||
write_expected_s2 expected_s2;
|
||||
write_psi_coef ~n_det:n_det psi_coef ;
|
||||
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
|
||||
write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
|
||||
;;
|
||||
|
||||
|
|
@ -298,6 +305,7 @@ Determinants ::
|
|||
n_int = %s
|
||||
bit_kind = %s
|
||||
n_det = %s
|
||||
n_states = %s
|
||||
expected_s2 = %s
|
||||
psi_coef = %s
|
||||
psi_det = %s
|
||||
|
|
@ -305,6 +313,7 @@ psi_det = %s
|
|||
(b.n_int |> N_int_number.to_string)
|
||||
(b.bit_kind |> Bit_kind.to_string)
|
||||
(b.n_det |> Det_number.to_string)
|
||||
(b.n_states |> States_number.to_string)
|
||||
(b.expected_s2 |> Positive_float.to_string)
|
||||
(b.psi_coef |> Array.to_list |> List.map ~f:Det_coef.to_string
|
||||
|> String.concat ~sep:", ")
|
||||
|
|
@ -433,14 +442,83 @@ psi_det = %s
|
|||
|> Bit_kind.to_int)
|
||||
and n_int =
|
||||
Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
|
||||
and n_states =
|
||||
Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ())
|
||||
in
|
||||
let s =
|
||||
String.concat [ header ; bitkind ; n_int ; psi_coef ; psi_det]
|
||||
String.concat [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det]
|
||||
in
|
||||
|
||||
|
||||
|
||||
|
||||
Generic_input_of_rst.evaluate_sexp t_of_sexp s
|
||||
;;
|
||||
|
||||
let update_ndet n_det_new =
|
||||
Printf.printf "Reducing n_det to %d\n" (Det_number.to_int n_det_new);
|
||||
let n_det_new =
|
||||
Det_number.to_int n_det_new
|
||||
in
|
||||
let det =
|
||||
read ()
|
||||
in
|
||||
let n_det_old, n_states =
|
||||
Det_number.to_int det.n_det,
|
||||
States_number.to_int det.n_states
|
||||
in
|
||||
if n_det_new = n_det_old then
|
||||
()
|
||||
;
|
||||
if n_det_new > n_det_new then
|
||||
failwith @@ Printf.sprintf "Requested n_det should be less than %d" n_det_old
|
||||
;
|
||||
for j=0 to (n_states-1) do
|
||||
let ishift_old, ishift_new =
|
||||
j*n_det_old,
|
||||
j*n_det_new
|
||||
in
|
||||
for i=0 to (n_det_new-1) do
|
||||
det.psi_coef.(i+ishift_new) <- det.psi_coef.(i+ishift_old)
|
||||
done
|
||||
done
|
||||
;
|
||||
let new_det =
|
||||
{ det with n_det = (Det_number.of_int n_det_new) }
|
||||
in
|
||||
write new_det
|
||||
;;
|
||||
|
||||
let extract_state istate =
|
||||
Printf.printf "Extracting state %d\n" (States_number.to_int istate);
|
||||
let det =
|
||||
read ()
|
||||
in
|
||||
let n_det, n_states =
|
||||
Det_number.to_int det.n_det,
|
||||
States_number.to_int det.n_states
|
||||
in
|
||||
if (States_number.to_int istate) > n_states then
|
||||
failwith "State to extract should not be greater than n_states"
|
||||
;
|
||||
let j =
|
||||
(States_number.to_int istate) - 1
|
||||
in
|
||||
begin
|
||||
if (j>0) then
|
||||
let ishift =
|
||||
j*n_det
|
||||
in
|
||||
for i=0 to (n_det-1) do
|
||||
det.psi_coef.(i) <- det.psi_coef.(i+ishift)
|
||||
done
|
||||
end;
|
||||
let new_det =
|
||||
{ det with n_states = (States_number.of_int 1) }
|
||||
in
|
||||
write new_det
|
||||
;;
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -85,7 +85,7 @@ module Xyz = struct
|
|||
let of_string s =
|
||||
let flush state accu number =
|
||||
let n =
|
||||
if (number = "") then 0
|
||||
if (number = "") then 1
|
||||
else (Int.of_string number)
|
||||
in
|
||||
match state with
|
||||
|
|
|
|||
|
|
@ -47,6 +47,14 @@ let debug str =
|
|||
let zmq_context =
|
||||
ZMQ.Context.create ()
|
||||
|
||||
let () =
|
||||
let nproc =
|
||||
match Sys.getenv "OMP_NUM_THREADS" with
|
||||
| Some m -> int_of_string m
|
||||
| None -> 2
|
||||
in
|
||||
ZMQ.Context.set_io_threads zmq_context nproc
|
||||
|
||||
|
||||
let bind_socket ~socket_type ~socket ~port =
|
||||
let rec loop = function
|
||||
|
|
|
|||
|
|
@ -21,6 +21,9 @@ let spec =
|
|||
~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)"
|
||||
+> anon ("(xyz_file|zmt_file)" %: file )
|
||||
|
||||
type element =
|
||||
| Element of Element.t
|
||||
| Int_elem of (Nucl_number.t * Element.t)
|
||||
|
||||
(** Handle dummy atoms placed on bonds *)
|
||||
let dummy_centers ~threshold ~molecule ~nuclei =
|
||||
|
|
@ -115,17 +118,14 @@ let run ?o b c d m p cart xyz_file =
|
|||
(* Open basis set channels *)
|
||||
let basis_channel element =
|
||||
let key =
|
||||
Element.to_string element
|
||||
match element with
|
||||
| Element e -> Element.to_string e
|
||||
| Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e)
|
||||
in
|
||||
match Hashtbl.find basis_table key with
|
||||
| Some in_channel ->
|
||||
in_channel
|
||||
| None ->
|
||||
let msg =
|
||||
Printf.sprintf "%s is not defined in basis %s.%!"
|
||||
(Element.to_long_string element) b ;
|
||||
in
|
||||
failwith msg
|
||||
| None -> raise Not_found
|
||||
in
|
||||
|
||||
let temp_filename =
|
||||
|
|
@ -189,12 +189,21 @@ let run ?o b c d m p cart xyz_file =
|
|||
| Some (key, basis) -> (*Aux basis *)
|
||||
begin
|
||||
let elem =
|
||||
Element.of_string key
|
||||
try
|
||||
Element (Element.of_string key)
|
||||
with Element.ElementError _ ->
|
||||
let result =
|
||||
match (String.split ~on:',' key) with
|
||||
| i :: k :: [] -> (Nucl_number.of_int @@ int_of_string i, Element.of_string k)
|
||||
| _ -> failwith "Expected format is int,Element:basis"
|
||||
in Int_elem result
|
||||
and basis =
|
||||
String.lowercase basis
|
||||
in
|
||||
let key =
|
||||
Element.to_string elem
|
||||
match elem with
|
||||
| Element e -> Element.to_string e
|
||||
| Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e)
|
||||
in
|
||||
let new_channel =
|
||||
fetch_channel basis
|
||||
|
|
@ -202,7 +211,13 @@ let run ?o b c d m p cart xyz_file =
|
|||
begin
|
||||
match Hashtbl.add basis_table ~key:key ~data:new_channel with
|
||||
| `Ok -> ()
|
||||
| `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem))
|
||||
| `Duplicate ->
|
||||
let e =
|
||||
match elem with
|
||||
| Element e -> e
|
||||
| Int_elem (_,e) -> e
|
||||
in
|
||||
failwith ("Duplicate definition of basis for "^(Element.to_long_string e))
|
||||
end
|
||||
end
|
||||
end;
|
||||
|
|
@ -537,7 +552,20 @@ let run ?o b c d m p cart xyz_file =
|
|||
| Element.X -> Element.H
|
||||
| e -> e
|
||||
in
|
||||
Basis.read_element (basis_channel x.Atom.element) i e
|
||||
let key =
|
||||
Int_elem (i,x.Atom.element)
|
||||
in
|
||||
try
|
||||
Basis.read_element (basis_channel key) i e
|
||||
with Not_found ->
|
||||
let key =
|
||||
Element x.Atom.element
|
||||
in
|
||||
try
|
||||
Basis.read_element (basis_channel key) i e
|
||||
with Not_found ->
|
||||
failwith (Printf.sprintf "Basis not found for atom %d (%s)" (Nucl_number.to_int i)
|
||||
(Element.to_string x.Atom.element) )
|
||||
with
|
||||
| End_of_file -> failwith
|
||||
("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
|
||||
|
|
@ -647,6 +675,7 @@ atoms are taken from the same basis set, otherwise specific elements can be
|
|||
defined as follows:
|
||||
|
||||
-b \"cc-pcvdz | H:cc-pvdz | C:6-31g\"
|
||||
-b \"cc-pvtz | 1,H:sto-3g | 3,H:6-31g\"
|
||||
|
||||
If a file with the same name as the basis set exists, this file will be read.
|
||||
Otherwise, the basis set is obtained from the database.
|
||||
|
|
|
|||
|
|
@ -41,8 +41,8 @@ subroutine run_selection_slave(thread,iproc,energy)
|
|||
if (done) then
|
||||
ctask = ctask - 1
|
||||
else
|
||||
integer :: i_generator, i_generator_start, i_generator_max, step, N
|
||||
read (task,*) i_generator_start, i_generator_max, step, N
|
||||
integer :: i_generator, N
|
||||
read (task,*) i_generator, N
|
||||
if(buf%N == 0) then
|
||||
! Only first time
|
||||
call create_selection_buffer(N, N*2, buf)
|
||||
|
|
@ -50,9 +50,7 @@ subroutine run_selection_slave(thread,iproc,energy)
|
|||
else
|
||||
if(N /= buf%N) stop "N changed... wtf man??"
|
||||
end if
|
||||
do i_generator=i_generator_start,i_generator_max,step
|
||||
call select_connected(i_generator,energy,pt2,buf)
|
||||
enddo
|
||||
call select_connected(i_generator,energy,pt2,buf)
|
||||
endif
|
||||
|
||||
if(done .or. ctask == size(task_id)) then
|
||||
|
|
|
|||
|
|
@ -671,13 +671,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
|
|||
if(mat(1, p1, p2) == 0d0) cycle
|
||||
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
|
||||
logical, external :: is_in_wavefunction
|
||||
if (is_in_wavefunction(det,N_int)) then
|
||||
stop 'is_in_wf'
|
||||
cycle
|
||||
endif
|
||||
|
||||
if (do_ddci) then
|
||||
integer, external :: is_a_two_holes_two_particles
|
||||
logical, external :: is_a_two_holes_two_particles
|
||||
if (is_a_two_holes_two_particles(det)) then
|
||||
cycle
|
||||
endif
|
||||
|
|
@ -1219,35 +1215,42 @@ subroutine ZMQ_selection(N_in, pt2)
|
|||
|
||||
implicit none
|
||||
|
||||
character*(512) :: task
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
integer, intent(in) :: N_in
|
||||
type(selection_buffer) :: b
|
||||
integer :: i, N
|
||||
integer, external :: omp_get_thread_num
|
||||
double precision, intent(out) :: pt2(N_states)
|
||||
integer, parameter :: maxtasks=10000
|
||||
|
||||
|
||||
N = max(N_in,1)
|
||||
if (.True.) then
|
||||
PROVIDE pt2_e0_denominator
|
||||
N = max(N_in,1)
|
||||
provide nproc
|
||||
call new_parallel_job(zmq_to_qp_run_socket,"selection")
|
||||
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
|
||||
call zmq_set_running(zmq_to_qp_run_socket)
|
||||
call create_selection_buffer(N, N*2, b)
|
||||
endif
|
||||
|
||||
integer :: i_generator, i_generator_start, i_generator_max, step
|
||||
character*(20*maxtasks) :: task
|
||||
task = ' '
|
||||
|
||||
step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num ))
|
||||
step = max(1,step)
|
||||
do i= 1, N_det_generators,step
|
||||
i_generator_start = i
|
||||
i_generator_max = min(i+step-1,N_det_generators)
|
||||
write(task,*) i_generator_start, i_generator_max, 1, N
|
||||
integer :: k
|
||||
k=0
|
||||
do i= 1, N_det_generators
|
||||
k = k+1
|
||||
write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
|
||||
k = k+20
|
||||
if (k>20*maxtasks) then
|
||||
k=0
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||
endif
|
||||
enddo
|
||||
if (k > 0) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||
end do
|
||||
endif
|
||||
call zmq_set_running(zmq_to_qp_run_socket)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
|
||||
i = omp_get_thread_num()
|
||||
|
|
@ -1264,6 +1267,7 @@ subroutine ZMQ_selection(N_in, pt2)
|
|||
if (s2_eig) then
|
||||
call make_s2_eigenfunction
|
||||
endif
|
||||
call save_wavefunction
|
||||
endif
|
||||
end subroutine
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
BEGIN_PROVIDER [integer, n_points_angular_grid]
|
||||
implicit none
|
||||
n_points_angular_grid = 50
|
||||
END_PROVIDER
|
||||
!BEGIN_PROVIDER [integer, n_points_integration_angular_lebedev]
|
||||
!implicit none
|
||||
!n_points_integration_angular_lebedev = 50
|
||||
!END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, n_points_radial_grid]
|
||||
implicit none
|
||||
|
|
@ -9,15 +9,15 @@ BEGIN_PROVIDER [integer, n_points_radial_grid]
|
|||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_angular_grid,3) ]
|
||||
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_angular_grid)]
|
||||
BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular_lebedev,3) ]
|
||||
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular_lebedev)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! weights and grid points for the integration on the angular variables on
|
||||
! the unit sphere centered on (0,0,0)
|
||||
! According to the LEBEDEV scheme
|
||||
END_DOC
|
||||
call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points)
|
||||
call cal_quad(n_points_integration_angular_lebedev, angular_quadrature_points,weights_angular_points)
|
||||
include 'constants.include.F'
|
||||
integer :: i
|
||||
double precision :: accu
|
||||
|
|
@ -63,7 +63,7 @@ END_PROVIDER
|
|||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid,n_points_radial_grid,nucl_num)]
|
||||
BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num)]
|
||||
BEGIN_DOC
|
||||
! points for integration over space
|
||||
END_DOC
|
||||
|
|
@ -79,7 +79,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid
|
|||
double precision :: x,r
|
||||
x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1]
|
||||
r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration
|
||||
do k = 1, n_points_angular_grid ! explicit values of the grid points centered around each atom
|
||||
do k = 1, n_points_integration_angular_lebedev ! explicit values of the grid points centered around each atom
|
||||
grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r
|
||||
grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * r
|
||||
grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r
|
||||
|
|
@ -88,7 +88,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid
|
|||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_DOC
|
||||
! Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
|
||||
! the "n" discrete variable represents the nucleis which in this array is represented by the last dimension
|
||||
|
|
@ -102,7 +102,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
|
|||
! run over all points in space
|
||||
do j = 1, nucl_num ! that are referred to each atom
|
||||
do k = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
|
||||
do l = 1, n_points_angular_grid ! for each angular point attached to the "jth" atom
|
||||
do l = 1, n_points_integration_angular_lebedev ! for each angular point attached to the "jth" atom
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
r(2) = grid_points_per_atom(2,l,k,j)
|
||||
r(3) = grid_points_per_atom(3,l,k,j)
|
||||
|
|
@ -123,8 +123,8 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
|
|||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
implicit none
|
||||
integer :: i,j,k,l,m
|
||||
double precision :: contrib
|
||||
|
|
@ -132,7 +132,7 @@ END_PROVIDER
|
|||
double precision :: aos_array(ao_num),mos_array(mo_tot_num)
|
||||
do j = 1, nucl_num
|
||||
do k = 1, n_points_radial_grid -1
|
||||
do l = 1, n_points_angular_grid
|
||||
do l = 1, n_points_integration_angular_lebedev
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
double precision :: accu
|
||||
integer :: i,j,k,l
|
||||
double precision :: x
|
||||
double precision :: integrand(n_points_angular_grid), weights(n_points_angular_grid)
|
||||
double precision :: integrand(n_points_integration_angular_lebedev), weights(n_points_integration_angular_lebedev)
|
||||
double precision :: f_average_angular_alpha,f_average_angular_beta
|
||||
double precision :: derivative_knowles_function,knowles_function
|
||||
|
||||
|
|
@ -12,7 +12,7 @@
|
|||
! according ot equation (6) of the paper of Becke (JCP, (88), 1988)
|
||||
! Here the m index is referred to the w_m(r) weight functions of equation (22)
|
||||
! Run over all points of integrations : there are
|
||||
! n_points_radial_grid (i) * n_points_angular_grid (k)
|
||||
! n_points_radial_grid (i) * n_points_integration_angular_lebedev (k)
|
||||
do j = 1, nucl_num
|
||||
integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
|
||||
integral_density_beta_knowles_becke_per_atom(j) = 0.d0
|
||||
|
|
@ -20,7 +20,7 @@
|
|||
! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
|
||||
f_average_angular_alpha = 0.d0
|
||||
f_average_angular_beta = 0.d0
|
||||
do k = 1, n_points_angular_grid
|
||||
do k = 1, n_points_integration_angular_lebedev
|
||||
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
|
||||
f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
|
||||
enddo
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ]
|
|||
! E0 in the denominator of the PT2
|
||||
END_DOC
|
||||
if (initialize_pt2_E0_denominator) then
|
||||
pt2_E0_denominator(1:N_states) = CI_electronic_energy(1:N_states)
|
||||
pt2_E0_denominator(1:N_states) = psi_energy(1:N_states)
|
||||
! pt2_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
|
||||
! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
|
||||
call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator')
|
||||
|
|
|
|||
|
|
@ -1,8 +1,7 @@
|
|||
program pt2_stoch
|
||||
implicit none
|
||||
initialize_pt2_E0_denominator = .False.
|
||||
read_wf = .True.
|
||||
SOFT_TOUCH initialize_pt2_E0_denominator read_wf
|
||||
SOFT_TOUCH read_wf
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
call run
|
||||
end
|
||||
|
|
@ -17,31 +16,23 @@ subroutine run
|
|||
integer :: n_det_before, to_select
|
||||
double precision :: threshold_davidson_in
|
||||
|
||||
double precision :: E_CI_before(N_states), relative_error
|
||||
double precision :: E_CI_before, relative_error
|
||||
|
||||
if (.true.) then
|
||||
call ezfio_get_full_ci_zmq_energy(E_CI_before(1))
|
||||
pt2_e0_denominator(:) = E_CI_before(1) - nuclear_repulsion
|
||||
SOFT_TOUCH pt2_e0_denominator read_wf
|
||||
endif
|
||||
allocate (pt2(N_states))
|
||||
pt2 = 0.d0
|
||||
|
||||
E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
relative_error = 1.d-6
|
||||
relative_error = 1.d-3
|
||||
call ZMQ_pt2(pt2, relative_error)
|
||||
print *, 'Final step'
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
do k=1,N_states
|
||||
print *, 'State', k
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E = ', E_CI_before
|
||||
print *, 'E+PT2 = ', E_CI_before+pt2
|
||||
print *, '-----'
|
||||
enddo
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E = ', E_CI_before
|
||||
print *, 'E+PT2 = ', E_CI_before+pt2
|
||||
print *, '-----'
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before+pt2(1))
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
|
|||
double precision, allocatable :: pt2_detail(:,:), comb(:)
|
||||
logical, allocatable :: computed(:)
|
||||
integer, allocatable :: tbc(:)
|
||||
integer :: i, j, Ncomb, generator_per_task, i_generator_end
|
||||
integer :: i, j, k, Ncomb, generator_per_task, i_generator_end
|
||||
integer, external :: pt2_find
|
||||
|
||||
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
|
|
@ -32,7 +32,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
|
|||
sum2above = 0d0
|
||||
Nabove = 0d0
|
||||
|
||||
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral
|
||||
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight
|
||||
|
||||
!call random_seed()
|
||||
|
||||
|
|
@ -69,18 +69,18 @@ subroutine ZMQ_pt2(pt2,relative_error)
|
|||
|
||||
do i=1,tbc(0)
|
||||
if(tbc(i) > fragment_first) then
|
||||
write(task(ipos:ipos+20),'(I9,X,I9,''|'')') 0, tbc(i)
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, tbc(i)
|
||||
ipos += 20
|
||||
if (ipos > 64000) then
|
||||
if (ipos > 63980) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
|
||||
ipos=1
|
||||
tasks = .True.
|
||||
endif
|
||||
else
|
||||
do j=1,fragment_count
|
||||
write(task(ipos:ipos+20),'(I9,X,I9,''|'')') j, tbc(i)
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, tbc(i)
|
||||
ipos += 20
|
||||
if (ipos > 64000) then
|
||||
if (ipos > 63980) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
|
||||
ipos=1
|
||||
tasks = .True.
|
||||
|
|
@ -108,7 +108,12 @@ subroutine ZMQ_pt2(pt2,relative_error)
|
|||
call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
|
||||
|
||||
else
|
||||
pt2(1) = sum(pt2_detail(1,:))
|
||||
pt2 = 0.d0
|
||||
do i=1,N_det_generators
|
||||
do k=1,N_states
|
||||
pt2(k) = pt2(k) + pt2_detail(k,i)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
tbc(0) = 0
|
||||
|
|
@ -117,6 +122,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
|
|||
endif
|
||||
end do
|
||||
|
||||
deallocate(pt2_detail, comb, computed, tbc)
|
||||
|
||||
end subroutine
|
||||
|
||||
|
|
@ -196,11 +202,15 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
|
|||
|
||||
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
|
||||
pt2_mwen(N_states, N_det_generators) )
|
||||
actually_computed(:) = computed(:)
|
||||
do i=1,N_det_generators
|
||||
actually_computed(i) = computed(i)
|
||||
enddo
|
||||
|
||||
parts_to_get(:) = 1
|
||||
if(fragment_first > 0) then
|
||||
parts_to_get(1:fragment_first) = fragment_count
|
||||
do i=1,fragment_first
|
||||
parts_to_get(i) = fragment_count
|
||||
enddo
|
||||
endif
|
||||
|
||||
do i=1,tbc(0)
|
||||
|
|
@ -223,7 +233,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
|
|||
pullLoop : do while (more == 1)
|
||||
call pull_pt2_results(zmq_socket_pull, Nindex, index, pt2_mwen, task_id, ntask)
|
||||
do i=1,Nindex
|
||||
pt2_detail(:, index(i)) += pt2_mwen(:,i)
|
||||
pt2_detail(1:N_states, index(i)) += pt2_mwen(1:N_states,i)
|
||||
parts_to_get(index(i)) -= 1
|
||||
if(parts_to_get(index(i)) < 0) then
|
||||
print *, i, index(i), parts_to_get(index(i)), Nindex
|
||||
|
|
@ -273,12 +283,11 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
|
|||
if (dabs(eqt/avg) < relative_error) then
|
||||
pt2(1) = avg
|
||||
! exit pullLoop
|
||||
else
|
||||
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
|
||||
endif
|
||||
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
|
||||
end if
|
||||
end do pullLoop
|
||||
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
|
||||
|
||||
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_pull_socket(zmq_socket_pull)
|
||||
|
|
@ -375,9 +384,9 @@ END_PROVIDER
|
|||
|
||||
subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
|
||||
implicit none
|
||||
integer, intent(inout) :: Ncomb
|
||||
double precision, intent(out) :: comb(Ncomb)
|
||||
integer, intent(inout) :: tbc(0:size_tbc)
|
||||
integer, intent(inout) :: Ncomb
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer :: i, j, last_full, dets(comb_teeth), tbc_save
|
||||
integer :: icount, n
|
||||
|
|
@ -515,12 +524,15 @@ end subroutine
|
|||
pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
|
||||
end do
|
||||
|
||||
pt2_weight = pt2_weight / pt2_cweight(N_det_generators)
|
||||
pt2_cweight = pt2_cweight / pt2_cweight(N_det_generators)
|
||||
do i=1,N_det_generators
|
||||
pt2_weight(i) = pt2_weight(i) / pt2_cweight(N_det_generators)
|
||||
pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(N_det_generators)
|
||||
enddo
|
||||
|
||||
norm_left = 1d0
|
||||
|
||||
comb_step = 1d0/dfloat(comb_teeth)
|
||||
first_det_of_comb = 1
|
||||
do i=1,N_det_generators
|
||||
if(pt2_weight(i)/norm_left < comb_step*.5d0) then
|
||||
first_det_of_comb = i
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
|
|||
integer :: index
|
||||
integer :: Nindex
|
||||
|
||||
allocate(pt2_detail(N_states, N_det))
|
||||
allocate(pt2_detail(N_states, N_det_generators))
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
zmq_socket_push = new_zmq_push_socket(thread)
|
||||
call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
|
||||
|
|
@ -101,7 +101,7 @@ subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntas
|
|||
implicit none
|
||||
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_push
|
||||
double precision, intent(in) :: pt2_detail(N_states, N_det)
|
||||
double precision, intent(in) :: pt2_detail(N_states, N_det_generators)
|
||||
integer, intent(in) :: ntask, N, index, task_id(*)
|
||||
integer :: rc
|
||||
|
||||
|
|
@ -133,7 +133,7 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
|
|||
use selection_types
|
||||
implicit none
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||
double precision, intent(inout) :: pt2_detail(N_states, N_det)
|
||||
double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
|
||||
integer, intent(out) :: index
|
||||
integer, intent(out) :: N, ntask, task_id(*)
|
||||
integer :: rc, rn, i
|
||||
|
|
@ -150,18 +150,22 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
|
|||
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
|
||||
if(rc /= 4) stop "pull"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id, ntask*4, 0)
|
||||
if(rc /= 4*ntask) stop "pull"
|
||||
|
||||
! Activate is zmq_socket_pull is a REP
|
||||
rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
|
||||
|
||||
do i=N+1,N_det_generators
|
||||
pt2_detail(1:N_states,i) = 0.d0
|
||||
enddo
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, pt2_workload, (N_det) ]
|
||||
BEGIN_PROVIDER [ double precision, pt2_workload, (N_det_generators) ]
|
||||
integer :: i
|
||||
do i=1,N_det
|
||||
pt2_workload(:) = dfloat(N_det - i + 1)**2
|
||||
do i=1,N_det_generators
|
||||
pt2_workload(i) = dfloat(N_det_generators - i + 1)**2
|
||||
end do
|
||||
pt2_workload = pt2_workload / sum(pt2_workload)
|
||||
END_PROVIDER
|
||||
|
|
|
|||
|
|
@ -26,7 +26,6 @@ subroutine run_selection_slave(thread,iproc,energy)
|
|||
call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
|
||||
if(worker_id == -1) then
|
||||
print *, "WORKER -1"
|
||||
!call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
return
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -41,31 +41,33 @@ subroutine sort_selection_buffer(b)
|
|||
implicit none
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
double precision, allocatable :: vals(:), absval(:)
|
||||
double precision, allocatable:: absval(:)
|
||||
integer, allocatable :: iorder(:)
|
||||
integer(bit_kind), allocatable :: detmp(:,:,:)
|
||||
double precision, pointer :: vals(:)
|
||||
integer(bit_kind), pointer :: detmp(:,:,:)
|
||||
integer :: i, nmwen
|
||||
logical, external :: detEq
|
||||
nmwen = min(b%N, b%cur)
|
||||
|
||||
|
||||
allocate(iorder(b%cur), detmp(N_int, 2, nmwen), absval(b%cur), vals(nmwen))
|
||||
allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
|
||||
absval = -dabs(b%val(:b%cur))
|
||||
do i=1,b%cur
|
||||
iorder(i) = i
|
||||
end do
|
||||
call dsort(absval, iorder, b%cur)
|
||||
|
||||
! Optimal for almost sorted data
|
||||
call insertion_dsort(absval, iorder, b%cur)
|
||||
do i=1, nmwen
|
||||
detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
|
||||
detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
|
||||
vals(i) = b%val(iorder(i))
|
||||
end do
|
||||
b%det = 0_bit_kind
|
||||
b%val = 0d0
|
||||
b%det(1:N_int,1,1:nmwen) = detmp(1:N_int,1,1:nmwen)
|
||||
b%det(1:N_int,2,1:nmwen) = detmp(1:N_int,2,1:nmwen)
|
||||
b%val(1:nmwen) = vals(1:nmwen)
|
||||
do i=nmwen+1, size(vals)
|
||||
vals(i) = 0.d0
|
||||
enddo
|
||||
deallocate(b%det, b%val)
|
||||
b%det => detmp
|
||||
b%val => vals
|
||||
b%mini = max(b%mini,dabs(b%val(b%N)))
|
||||
b%cur = nmwen
|
||||
end subroutine
|
||||
|
|
|
|||
|
|
@ -1,9 +1,9 @@
|
|||
module selection_types
|
||||
type selection_buffer
|
||||
integer :: N, cur
|
||||
integer(8), allocatable :: det(:,:,:)
|
||||
double precision, allocatable :: val(:)
|
||||
double precision :: mini
|
||||
integer(8) , pointer :: det(:,:,:)
|
||||
double precision, pointer :: val(:)
|
||||
double precision :: mini
|
||||
endtype
|
||||
end module
|
||||
|
||||
|
|
|
|||
|
|
@ -10,26 +10,38 @@ subroutine ZMQ_selection(N_in, pt2)
|
|||
integer :: i, N
|
||||
integer, external :: omp_get_thread_num
|
||||
double precision, intent(out) :: pt2(N_states)
|
||||
integer, parameter :: maxtasks=10000
|
||||
|
||||
|
||||
PROVIDE fragment_count
|
||||
|
||||
N = max(N_in,1)
|
||||
if (.True.) then
|
||||
PROVIDE pt2_e0_denominator
|
||||
N = max(N_in,1)
|
||||
provide nproc
|
||||
call new_parallel_job(zmq_to_qp_run_socket,"selection")
|
||||
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
|
||||
call zmq_set_running(zmq_to_qp_run_socket)
|
||||
call create_selection_buffer(N, N*2, b)
|
||||
endif
|
||||
|
||||
character(len=:), allocatable :: task
|
||||
task = repeat(' ',20*N_det_generators)
|
||||
character*(20*maxtasks) :: task
|
||||
task = ' '
|
||||
|
||||
integer :: k
|
||||
k=0
|
||||
do i= 1, N_det_generators
|
||||
write(task(20*(i-1)+1:20*i),'(I9,X,I9,''|'')') i, N
|
||||
k = k+1
|
||||
write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
|
||||
k = k+20
|
||||
if (k>20*maxtasks) then
|
||||
k=0
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||
endif
|
||||
end do
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||
if (k > 0) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||
endif
|
||||
call zmq_set_running(zmq_to_qp_run_socket)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
|
||||
i = omp_get_thread_num()
|
||||
|
|
@ -48,6 +60,7 @@ subroutine ZMQ_selection(N_in, pt2)
|
|||
endif
|
||||
call save_wavefunction
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
|
|
@ -83,7 +96,7 @@ subroutine selection_collector(b, pt2)
|
|||
real :: time, time0
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
zmq_socket_pull = new_zmq_pull_socket()
|
||||
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
|
||||
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators))
|
||||
done = 0
|
||||
more = 1
|
||||
pt2(:) = 0d0
|
||||
|
|
|
|||
|
|
@ -23,33 +23,39 @@
|
|||
allocate(pathTo(N_det_non_ref))
|
||||
|
||||
pathTo(:) = 0
|
||||
is_active_exc(:) = .false.
|
||||
is_active_exc(:) = .True.
|
||||
n_exc_active = 0
|
||||
|
||||
do hh = 1, hh_shortcut(0)
|
||||
do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
|
||||
do II = 1, N_det_ref
|
||||
! do hh = 1, hh_shortcut(0)
|
||||
! do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
|
||||
! do II = 1, N_det_ref
|
||||
!
|
||||
! call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
|
||||
! if(.not. ok) cycle
|
||||
!
|
||||
! call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
|
||||
! if(.not. ok) cycle
|
||||
!
|
||||
! ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
|
||||
! if(ind == -1) cycle
|
||||
!
|
||||
! logical, external :: is_a_two_holes_two_particles
|
||||
! if (is_a_two_holes_two_particles(myDet)) then
|
||||
! is_active_exc(pp) = .False.
|
||||
! endif
|
||||
|
||||
call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
|
||||
if(.not. ok) cycle
|
||||
! ind = psi_non_ref_sorted_idx(ind)
|
||||
! if(pathTo(ind) == 0) then
|
||||
! pathTo(ind) = pp
|
||||
! else
|
||||
! is_active_exc(pp) = .true.
|
||||
! is_active_exc(pathTo(ind)) = .true.
|
||||
! end if
|
||||
|
||||
call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
|
||||
if(.not. ok) cycle
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
|
||||
ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
|
||||
if(ind == -1) cycle
|
||||
|
||||
ind = psi_non_ref_sorted_idx(ind)
|
||||
if(pathTo(ind) == 0) then
|
||||
pathTo(ind) = pp
|
||||
else
|
||||
is_active_exc(pp) = .true.
|
||||
is_active_exc(pathTo(ind)) = .true.
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!is_active_exc=.true.
|
||||
do hh = 1, hh_shortcut(0)
|
||||
do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
|
||||
if(is_active_exc(pp)) then
|
||||
|
|
@ -66,6 +72,32 @@
|
|||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ logical, has_a_unique_parent, (N_det_non_ref) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! True if the determinant in the non-reference has a unique parent
|
||||
END_DOC
|
||||
integer :: i,j,n
|
||||
integer :: degree
|
||||
do j=1,N_det_non_ref
|
||||
has_a_unique_parent(j) = .True.
|
||||
n=0
|
||||
do i=1,N_det_ref
|
||||
call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,j), degree, N_int)
|
||||
if (degree < 2) then
|
||||
n = n+1
|
||||
if (n > 1) then
|
||||
has_a_unique_parent(j) = .False.
|
||||
exit
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, n_exc_active_sze ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
|
@ -96,7 +128,7 @@ END_PROVIDER
|
|||
!$OMP active_excitation_to_determinants_val, active_excitation_to_determinants_idx)&
|
||||
!$OMP shared(hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, &
|
||||
!$OMP psi_non_ref_sorted_idx, psi_ref, N_det_ref, N_states)&
|
||||
!$OMP shared(is_active_exc, active_hh_idx, active_pp_idx, n_exc_active)&
|
||||
!$OMP shared(active_hh_idx, active_pp_idx, n_exc_active)&
|
||||
!$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh, s)
|
||||
allocate(lref(N_det_non_ref))
|
||||
!$OMP DO schedule(dynamic)
|
||||
|
|
|
|||
|
|
@ -351,11 +351,11 @@ logical function is_generable(det1, det2, Nint)
|
|||
integer, intent(in) :: Nint
|
||||
integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
|
||||
integer :: degree, f, exc(0:2, 2, 2), t
|
||||
integer*2 :: h1, h2, p1, p2, s1, s2
|
||||
integer :: h1, h2, p1, p2, s1, s2
|
||||
integer, external :: searchExc
|
||||
logical, external :: excEq
|
||||
double precision :: phase
|
||||
integer*2 :: tmp_array(4)
|
||||
integer :: tmp_array(4)
|
||||
|
||||
is_generable = .false.
|
||||
call get_excitation(det1, det2, exc, degree, phase, Nint)
|
||||
|
|
@ -366,7 +366,7 @@ logical function is_generable(det1, det2, Nint)
|
|||
end if
|
||||
if(degree > 2) stop "?22??"
|
||||
|
||||
call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
|
||||
if(degree == 1) then
|
||||
h2 = h1
|
||||
|
|
@ -454,7 +454,7 @@ integer function searchExc(excs, exc, n)
|
|||
use bitmasks
|
||||
|
||||
integer, intent(in) :: n
|
||||
integer*2,intent(in) :: excs(4,n), exc(4)
|
||||
integer,intent(in) :: excs(4,n), exc(4)
|
||||
integer :: l, h, c
|
||||
integer, external :: excCmp
|
||||
logical, external :: excEq
|
||||
|
|
@ -519,8 +519,8 @@ subroutine sort_exc(key, N_key)
|
|||
|
||||
|
||||
integer, intent(in) :: N_key
|
||||
integer*2,intent(inout) :: key(4,N_key)
|
||||
integer*2 :: tmp(4)
|
||||
integer,intent(inout) :: key(4,N_key)
|
||||
integer :: tmp(4)
|
||||
integer :: i,ni
|
||||
|
||||
|
||||
|
|
@ -542,7 +542,7 @@ end subroutine
|
|||
|
||||
logical function exc_inf(exc1, exc2)
|
||||
implicit none
|
||||
integer*2,intent(in) :: exc1(4), exc2(4)
|
||||
integer,intent(in) :: exc1(4), exc2(4)
|
||||
integer :: i
|
||||
exc_inf = .false.
|
||||
do i=1,4
|
||||
|
|
@ -564,9 +564,9 @@ subroutine tamise_exc(key, no, n, N_key)
|
|||
! Uncodumented : TODO
|
||||
END_DOC
|
||||
integer,intent(in) :: no, n, N_key
|
||||
integer*2,intent(inout) :: key(4, N_key)
|
||||
integer,intent(inout) :: key(4, N_key)
|
||||
integer :: k,j
|
||||
integer*2 :: tmp(4)
|
||||
integer :: tmp(4)
|
||||
logical :: exc_inf
|
||||
integer :: ni
|
||||
|
||||
|
|
@ -595,8 +595,9 @@ end subroutine
|
|||
|
||||
subroutine dec_exc(exc, h1, h2, p1, p2)
|
||||
implicit none
|
||||
integer :: exc(0:2,2,2), s1, s2, degree
|
||||
integer*2, intent(out) :: h1, h2, p1, p2
|
||||
integer, intent(in) :: exc(0:2,2,2)
|
||||
integer, intent(out) :: h1, h2, p1, p2
|
||||
integer :: degree, s1, s2
|
||||
|
||||
degree = exc(0,1,1) + exc(0,1,2)
|
||||
|
||||
|
|
@ -607,7 +608,7 @@ subroutine dec_exc(exc, h1, h2, p1, p2)
|
|||
|
||||
if(degree == 0) return
|
||||
|
||||
call decode_exc_int2(exc, degree, h1, p1, h2, p2, s1, s2)
|
||||
call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
|
||||
|
||||
h1 += mo_tot_num * (s1-1)
|
||||
p1 += mo_tot_num * (s1-1)
|
||||
|
|
@ -639,7 +640,7 @@ end subroutine
|
|||
&BEGIN_PROVIDER [ integer, N_ex_exists ]
|
||||
implicit none
|
||||
integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
|
||||
integer*2 :: h1, h2, p1, p2
|
||||
integer :: h1, h2, p1, p2
|
||||
double precision :: phase
|
||||
logical,allocatable :: hh(:,:) , pp(:,:)
|
||||
|
||||
|
|
@ -739,12 +740,12 @@ END_PROVIDER
|
|||
double precision :: phase
|
||||
|
||||
|
||||
double precision, allocatable :: rho_mrcc_init(:)
|
||||
double precision, allocatable :: rho_mrcc_inact(:)
|
||||
integer :: a_coll, at_roww
|
||||
|
||||
print *, "TI", hh_nex, N_det_non_ref
|
||||
|
||||
allocate(rho_mrcc_init(N_det_non_ref))
|
||||
allocate(rho_mrcc_inact(N_det_non_ref))
|
||||
allocate(x_new(hh_nex))
|
||||
allocate(x(hh_nex), AtB(hh_nex))
|
||||
|
||||
|
|
@ -756,7 +757,7 @@ END_PROVIDER
|
|||
!$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
|
||||
!$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx)
|
||||
|
||||
!$OMP DO schedule(dynamic, 100)
|
||||
!$OMP DO schedule(static, 100)
|
||||
do at_roww = 1, n_exc_active ! hh_nex
|
||||
at_row = active_pp_idx(at_roww)
|
||||
do i=1,active_excitation_to_determinants_idx(0,at_roww)
|
||||
|
|
@ -775,7 +776,7 @@ END_PROVIDER
|
|||
X(a_col) = AtB(a_col)
|
||||
end do
|
||||
|
||||
rho_mrcc_init = 0d0
|
||||
rho_mrcc_inact(:) = 0d0
|
||||
|
||||
allocate(lref(N_det_ref))
|
||||
do hh = 1, hh_shortcut(0)
|
||||
|
|
@ -799,19 +800,15 @@ END_PROVIDER
|
|||
X(pp) = AtB(pp)
|
||||
do II=1,N_det_ref
|
||||
if(lref(II) > 0) then
|
||||
rho_mrcc_init(lref(II)) = psi_ref_coef(II,s) * X(pp)
|
||||
rho_mrcc_inact(lref(II)) = psi_ref_coef(II,s) * X(pp)
|
||||
else if(lref(II) < 0) then
|
||||
rho_mrcc_init(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
|
||||
rho_mrcc_inact(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
deallocate(lref)
|
||||
|
||||
do i=1,N_det_non_ref
|
||||
rho_mrcc(i,s) = rho_mrcc_init(i)
|
||||
enddo
|
||||
|
||||
x_new = x
|
||||
|
||||
double precision :: factor, resold
|
||||
|
|
@ -839,7 +836,10 @@ END_PROVIDER
|
|||
print *, k, res, 1.d0 - res/resold
|
||||
endif
|
||||
|
||||
if ( (res < 1d-10).or.(res/resold > 0.99d0) ) then
|
||||
if ( res < 1d-10 ) then
|
||||
exit
|
||||
endif
|
||||
if ( (res/resold > 0.99d0) ) then
|
||||
exit
|
||||
endif
|
||||
resold = res
|
||||
|
|
@ -848,38 +848,60 @@ END_PROVIDER
|
|||
dIj_unique(1:size(X), s) = X(1:size(X))
|
||||
print *, k, res, 1.d0 - res/resold
|
||||
|
||||
enddo
|
||||
|
||||
do s=1,N_states
|
||||
do i=1,N_det_non_ref
|
||||
rho_mrcc(i,s) = 0.d0
|
||||
enddo
|
||||
|
||||
do a_coll=1,n_exc_active
|
||||
a_col = active_pp_idx(a_coll)
|
||||
do j=1,N_det_non_ref
|
||||
i = active_excitation_to_determinants_idx(j,a_coll)
|
||||
if (i==0) exit
|
||||
if (rho_mrcc_inact(i) /= 0.d0) then
|
||||
call debug_det(psi_non_ref(1,1,i),N_int)
|
||||
stop
|
||||
endif
|
||||
rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s)
|
||||
enddo
|
||||
end do
|
||||
|
||||
norm = 0.d0
|
||||
do i=1,N_det_non_ref
|
||||
norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
|
||||
enddo
|
||||
! Norm now contains the norm of A.X
|
||||
|
||||
double precision :: norm2_ref, norm2_inact, a, b, c, Delta
|
||||
! Psi = Psi_ref + Psi_inactive + f*Psi_active
|
||||
! Find f to normalize Psi
|
||||
|
||||
norm2_ref = 0.d0
|
||||
do i=1,N_det_ref
|
||||
norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
|
||||
norm2_ref = norm2_ref + psi_ref_coef(i,s)*psi_ref_coef(i,s)
|
||||
enddo
|
||||
! Norm now contains the norm of Psi + A.X
|
||||
|
||||
|
||||
a = 0.d0
|
||||
do i=1,N_det_non_ref
|
||||
a = a + rho_mrcc(i,s)*rho_mrcc(i,s)
|
||||
enddo
|
||||
|
||||
norm = a + norm2_ref
|
||||
print *, "norm : ", sqrt(norm)
|
||||
enddo
|
||||
|
||||
norm = sqrt((1.d0-norm2_ref)/a)
|
||||
|
||||
! Renormalize Psi+A.X
|
||||
do i=1,N_det_non_ref
|
||||
rho_mrcc(i,s) = rho_mrcc(i,s) * norm
|
||||
enddo
|
||||
|
||||
!norm = norm2_ref
|
||||
!do i=1,N_det_non_ref
|
||||
! norm = norm + rho_mrcc(i,s)**2
|
||||
!enddo
|
||||
!print *, 'check', norm
|
||||
!stop
|
||||
|
||||
|
||||
|
||||
do s=1,N_states
|
||||
norm = 0.d0
|
||||
double precision :: f, g, gmax
|
||||
gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s)))
|
||||
gmax = maxval(dabs(psi_non_ref_coef(:,s)))
|
||||
do i=1,N_det_non_ref
|
||||
if (lambda_type == 2) then
|
||||
f = 1.d0
|
||||
|
|
@ -891,41 +913,22 @@ END_PROVIDER
|
|||
f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
|
||||
|
||||
! Avoid numerical instabilities
|
||||
! g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) )
|
||||
g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax))
|
||||
f = min(f, g)
|
||||
f = max(f,-g)
|
||||
|
||||
endif
|
||||
|
||||
norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
|
||||
norm = norm + (rho_mrcc(i,s)*f)**2
|
||||
rho_mrcc(i,s) = f
|
||||
enddo
|
||||
! norm now contains the norm of |T.Psi_0>
|
||||
! rho_mrcc now contains the f factors
|
||||
|
||||
f = 1.d0/norm
|
||||
! f now contains 1/ <T.Psi_0|T.Psi_0>
|
||||
|
||||
norm = 0.d0
|
||||
do i=1,N_det_non_ref
|
||||
norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
|
||||
enddo
|
||||
! norm now contains <Psi_SD|Psi_SD>
|
||||
f = dsqrt(f*norm)
|
||||
! f normalises T.Psi_0 such that (1+T)|Psi> is normalized
|
||||
! rho_mrcc now contains the mu_i factors
|
||||
|
||||
print *, 'norm of |T Psi_0> = ', dsqrt(norm)
|
||||
norm = norm*f
|
||||
if (dsqrt(norm) > 1.d0) then
|
||||
if (norm > 1.d0) then
|
||||
stop 'Error : Norm of the SD larger than the norm of the reference.'
|
||||
endif
|
||||
|
||||
do i=1,N_det_non_ref
|
||||
rho_mrcc(i,s) = rho_mrcc(i,s) * f
|
||||
enddo
|
||||
! rho_mrcc now contains the product of the scaling factors and the
|
||||
! normalization constant
|
||||
|
||||
end do
|
||||
|
||||
END_PROVIDER
|
||||
|
|
@ -1028,11 +1031,11 @@ double precision function get_dij(det1, det2, s, Nint)
|
|||
integer, intent(in) :: s, Nint
|
||||
integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
|
||||
integer :: degree, f, exc(0:2, 2, 2), t
|
||||
integer*2 :: h1, h2, p1, p2, s1, s2
|
||||
integer :: h1, h2, p1, p2, s1, s2
|
||||
integer, external :: searchExc
|
||||
logical, external :: excEq
|
||||
double precision :: phase
|
||||
integer*2 :: tmp_array(4)
|
||||
integer :: tmp_array(4)
|
||||
|
||||
get_dij = 0d0
|
||||
call get_excitation(det1, det2, exc, degree, phase, Nint)
|
||||
|
|
@ -1041,7 +1044,7 @@ double precision function get_dij(det1, det2, s, Nint)
|
|||
stop "get_dij"
|
||||
end if
|
||||
|
||||
call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
|
||||
if(degree == 1) then
|
||||
h2 = h1
|
||||
|
|
@ -1074,8 +1077,8 @@ double precision function get_dij(det1, det2, s, Nint)
|
|||
end function
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer*2, hh_exists, (4, N_hh_exists) ]
|
||||
&BEGIN_PROVIDER [ integer*2, pp_exists, (4, N_pp_exists) ]
|
||||
BEGIN_PROVIDER [ integer, hh_exists, (4, N_hh_exists) ]
|
||||
&BEGIN_PROVIDER [ integer, pp_exists, (4, N_pp_exists) ]
|
||||
&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_hh_exists + 1) ]
|
||||
&BEGIN_PROVIDER [ integer, hh_nex ]
|
||||
implicit none
|
||||
|
|
@ -1090,9 +1093,9 @@ end function
|
|||
! hh_nex : Total number of excitation operators
|
||||
!
|
||||
END_DOC
|
||||
integer*2,allocatable :: num(:,:)
|
||||
integer,allocatable :: num(:,:)
|
||||
integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
|
||||
integer*2 :: h1, h2, p1, p2
|
||||
integer :: h1, h2, p1, p2
|
||||
double precision :: phase
|
||||
logical, external :: excEq
|
||||
|
||||
|
|
@ -1118,19 +1121,19 @@ end function
|
|||
|
||||
hh_shortcut(0) = 1
|
||||
hh_shortcut(1) = 1
|
||||
hh_exists(:,1) = (/1_2, num(1,1), 1_2, num(2,1)/)
|
||||
pp_exists(:,1) = (/1_2, num(3,1), 1_2, num(4,1)/)
|
||||
hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/)
|
||||
pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/)
|
||||
s = 1
|
||||
do i=2,n
|
||||
if(.not. excEq(num(1,i), num(1,s))) then
|
||||
s += 1
|
||||
num(:, s) = num(:, i)
|
||||
pp_exists(:,s) = (/1_2, num(3,s), 1_2, num(4,s)/)
|
||||
pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/)
|
||||
if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. &
|
||||
hh_exists(4, hh_shortcut(0)) /= num(2,s)) then
|
||||
hh_shortcut(0) += 1
|
||||
hh_shortcut(hh_shortcut(0)) = s
|
||||
hh_exists(:,hh_shortcut(0)) = (/1_2, num(1,s), 1_2, num(2,s)/)
|
||||
hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
|
@ -1178,7 +1181,7 @@ END_PROVIDER
|
|||
|
||||
logical function excEq(exc1, exc2)
|
||||
implicit none
|
||||
integer*2, intent(in) :: exc1(4), exc2(4)
|
||||
integer, intent(in) :: exc1(4), exc2(4)
|
||||
integer :: i
|
||||
excEq = .false.
|
||||
do i=1, 4
|
||||
|
|
@ -1190,7 +1193,7 @@ end function
|
|||
|
||||
integer function excCmp(exc1, exc2)
|
||||
implicit none
|
||||
integer*2, intent(in) :: exc1(4), exc2(4)
|
||||
integer, intent(in) :: exc1(4), exc2(4)
|
||||
integer :: i
|
||||
excCmp = 0
|
||||
do i=1, 4
|
||||
|
|
@ -1209,8 +1212,8 @@ subroutine apply_hole_local(det, exc, res, ok, Nint)
|
|||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer*2, intent(in) :: exc(4)
|
||||
integer*2 :: s1, s2, h1, h2
|
||||
integer, intent(in) :: exc(4)
|
||||
integer :: s1, s2, h1, h2
|
||||
integer(bit_kind),intent(in) :: det(Nint, 2)
|
||||
integer(bit_kind),intent(out) :: res(Nint, 2)
|
||||
logical, intent(out) :: ok
|
||||
|
|
@ -1246,8 +1249,8 @@ subroutine apply_particle_local(det, exc, res, ok, Nint)
|
|||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer*2, intent(in) :: exc(4)
|
||||
integer*2 :: s1, s2, p1, p2
|
||||
integer, intent(in) :: exc(4)
|
||||
integer :: s1, s2, p1, p2
|
||||
integer(bit_kind),intent(in) :: det(Nint, 2)
|
||||
integer(bit_kind),intent(out) :: res(Nint, 2)
|
||||
logical, intent(out) :: ok
|
||||
|
|
|
|||
|
|
@ -1167,11 +1167,11 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
|
|||
print*, 'e corr perturb EN',accu(state_target)
|
||||
print*, ''
|
||||
print*, 'coef diagonalized'
|
||||
write(*,'(100(F16.10,X))')psi_in_out_coef(:,state_target)
|
||||
write(*,'(100(F16.10,1X))')psi_in_out_coef(:,state_target)
|
||||
print*, 'coef_perturb'
|
||||
write(*,'(100(F16.10,X))')coef_perturb(:)
|
||||
write(*,'(100(F16.10,1X))')coef_perturb(:)
|
||||
print*, 'coef_perturb EN'
|
||||
write(*,'(100(F16.10,X))')coef_perturb_bis(:)
|
||||
write(*,'(100(F16.10,1X))')coef_perturb_bis(:)
|
||||
endif
|
||||
integer :: k
|
||||
do k = 1, N_det_ref
|
||||
|
|
|
|||
|
|
@ -22,7 +22,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
|
|||
|
||||
integer :: elec_num_tab_local(2)
|
||||
integer :: i,j,accu_elec,k
|
||||
integer :: det_tmp(N_int), det_tmp_bis(N_int)
|
||||
integer(bit_kind) :: det_tmp(N_int), det_tmp_bis(N_int)
|
||||
double precision :: phase
|
||||
double precision :: norm_factor
|
||||
! print*, orb,hole_particle,spin_exc
|
||||
|
|
|
|||
|
|
@ -210,10 +210,6 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
|
|||
|
||||
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
|
||||
hab = (fock_operator_local(aorb,borb,kspin) ) * phase
|
||||
if(isnan(hab))then
|
||||
print*, '1'
|
||||
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)
|
||||
if(ispin == jspin)then
|
||||
|
|
@ -255,7 +251,8 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
|
|||
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,kspin) * phase
|
||||
if(isnan(hab))then
|
||||
! if(isnan(hab))then
|
||||
if(hab /= hab)then
|
||||
print*, '2'
|
||||
stop
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@
|
|||
z_min = 0.d0
|
||||
z_max = 10.d0
|
||||
delta_z = 0.005d0
|
||||
N_z_pts = (z_max - z_min)/delta_z
|
||||
N_z_pts = int( (z_max - z_min)/delta_z )
|
||||
print*,'N_z_pts = ',N_z_pts
|
||||
|
||||
END_PROVIDER
|
||||
|
|
|
|||
|
|
@ -151,7 +151,7 @@ subroutine print_hcc
|
|||
integer :: i,j
|
||||
print*,'Z AU GAUSS MHZ cm^-1'
|
||||
do i = 1, nucl_num
|
||||
write(*,'(I2,X,F4.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)
|
||||
write(*,'(I2,1X,F4.1,1X,4(F16.6,1X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
|
|
|||
|
|
@ -126,7 +126,7 @@ subroutine print_mulliken_sd
|
|||
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)
|
||||
write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,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
|
||||
|
|
@ -142,7 +142,7 @@ subroutine print_mulliken_sd
|
|||
accu = 0.d0
|
||||
do i = 0, ao_l_max
|
||||
accu += spin_population_angular_momentum_per_atom(i,j)
|
||||
write(*,'(XX,I3,XX,A4,X,A4,X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
|
||||
write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
|
||||
print*,'sum = ',accu
|
||||
enddo
|
||||
enddo
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ subroutine run
|
|||
write(*,'(A)') '============='
|
||||
write(*,'(A)') ''
|
||||
do istate=1,N_states
|
||||
call get_occupation_from_dets(occupation,1)
|
||||
call get_occupation_from_dets(occupation,istate)
|
||||
write(*,'(A)') ''
|
||||
write(*,'(A,I3)'), 'State ', istate
|
||||
write(*,'(A)') '---------------'
|
||||
|
|
@ -30,13 +30,13 @@ subroutine run
|
|||
if (occupation(i) > 1.999d0) then
|
||||
class(0,1) += 1
|
||||
class( class(0,1), 1) = i
|
||||
else if (occupation(i) > 1.95d0) then
|
||||
else if (occupation(i) > 1.97d0) then
|
||||
class(0,2) += 1
|
||||
class( class(0,2), 2) = i
|
||||
else if (occupation(i) < 0.001d0) then
|
||||
class(0,5) += 1
|
||||
class( class(0,5), 5) = i
|
||||
else if (occupation(i) < 0.01d0) then
|
||||
else if (occupation(i) < 0.03d0) then
|
||||
class(0,4) += 1
|
||||
class( class(0,4), 4) = i
|
||||
else
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
! DO NOT MODIFY BY HAND
|
||||
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
|
||||
! from file /ccc/work/cont003/gen1738/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
|
||||
! from file /panfs/panasas/cnt0024/cpq1738/scemama/workdir/quantum_package/src/mrcc_selected/EZFIO.cfg
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, thresh_dressed_ci ]
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ interface: ezfio
|
|||
type: logical
|
||||
doc: Compute perturbative contribution of the Triples
|
||||
interface: ezfio,provider,ocaml
|
||||
default: true
|
||||
default: false
|
||||
|
||||
[energy]
|
||||
type: double precision
|
||||
|
|
|
|||
|
|
@ -38,7 +38,7 @@ use bitmasks
|
|||
do p=hh_shortcut(h), hh_shortcut(h+1)-1
|
||||
call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int)
|
||||
if(ok) n = n + 1
|
||||
if(n > N_det_non_ref) stop "MRCC..."
|
||||
if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
|
||||
end do
|
||||
n = n - 1
|
||||
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
convert output of gamess/GAU$$IAN to ezfio
|
||||
|
||||
Usage:
|
||||
qp_convert_output_to_ezfio.py <file.out> [--ezfio=<folder.ezfio>]
|
||||
qp_convert_output_to_ezfio.py <file.out> [-o <ezfio_directory>]
|
||||
|
||||
Option:
|
||||
file.out is the file to check (like gamess.out)
|
||||
|
|
@ -272,7 +272,7 @@ def write_ezfio(res, filename):
|
|||
#
|
||||
|
||||
# INPUT
|
||||
# {% for lanel,zcore, l_block in l_atom $}
|
||||
# {% for label,zcore, l_block in l_atom $}
|
||||
# #local l_block l=0}
|
||||
# {label} GEN {zcore} {len(l_block)-1 #lmax_block}
|
||||
# {% for l_param in l_block%}
|
||||
|
|
@ -280,6 +280,7 @@ def write_ezfio(res, filename):
|
|||
# {% for coef,n,zeta for l_param}
|
||||
# {coef,n, zeta}
|
||||
|
||||
|
||||
# OUTPUT
|
||||
|
||||
# Local are 1 array padded by max(n_max_block) when l == 0 (output:k_loc_max)
|
||||
|
|
@ -309,8 +310,16 @@ def write_ezfio(res, filename):
|
|||
array_l_max_block.append(l_max_block)
|
||||
array_z_remove.append(z_remove)
|
||||
|
||||
matrix.append([[coef_n_zeta.split()[1:] for coef_n_zeta in l.split('\n')] for l in array_party[1:]])
|
||||
|
||||
x = [[coef_n_zeta.split() for coef_n_zeta in l.split('\n')] \
|
||||
for l in array_party[1:] ]
|
||||
x = []
|
||||
for l in array_party[1:]:
|
||||
y = []
|
||||
for coef_n_zeta in l.split('\n'):
|
||||
z = coef_n_zeta.split()
|
||||
if z : y.append(z)
|
||||
x.append(y)
|
||||
matrix.append(x)
|
||||
return (matrix, array_l_max_block, array_z_remove)
|
||||
|
||||
def get_local_stuff(matrix):
|
||||
|
|
@ -319,7 +328,6 @@ def write_ezfio(res, filename):
|
|||
k_loc_max = max(len(i) for i in matrix_local_unpad)
|
||||
|
||||
matrix_local = [ pad(ll, k_loc_max, [0., 2, 0.]) for ll in matrix_local_unpad]
|
||||
|
||||
m_coef = [[float(i[0]) for i in atom] for atom in matrix_local]
|
||||
m_n = [[int(i[1]) - 2 for i in atom] for atom in matrix_local]
|
||||
m_zeta = [[float(i[2]) for i in atom] for atom in matrix_local]
|
||||
|
|
@ -343,9 +351,20 @@ def write_ezfio(res, filename):
|
|||
return (l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc)
|
||||
|
||||
try:
|
||||
pseudo_str = res_file.get_pseudo()
|
||||
pseudo_str = []
|
||||
label = ezfio.get_nuclei_nucl_label()
|
||||
for ecp in res.pseudo:
|
||||
pseudo_str += [ "%(label)s GEN %(zcore)d %(lmax)d" % { "label": label[ ecp["atom"]-1 ],
|
||||
"zcore": ecp["zcore"], "lmax": ecp["lmax"] } ]
|
||||
lmax = ecp["lmax"]
|
||||
for l in [lmax] + list(range(0,lmax)):
|
||||
pseudo_str += [ "%d"%len(ecp[str(l)]) ]
|
||||
for t in ecp[str(l)]:
|
||||
pseudo_str += [ "%f %d %f"%t ]
|
||||
pseudo_str += [""]
|
||||
pseudo_str = "\n".join(pseudo_str)
|
||||
|
||||
matrix, array_l_max_block, array_z_remove = parse_str(pseudo_str)
|
||||
|
||||
except:
|
||||
ezfio.set_pseudo_do_pseudo(False)
|
||||
else:
|
||||
|
|
@ -359,10 +378,12 @@ def write_ezfio(res, filename):
|
|||
ezfio.nuclei_nucl_charge = [i - j for i, j in zip(ezfio.nuclei_nucl_charge, array_z_remove)]
|
||||
|
||||
import math
|
||||
num_elec = sum(ezfio.nuclei_nucl_charge)
|
||||
num_elec_diff = sum(array_z_remove)/2
|
||||
nalpha = ezfio.get_electrons_elec_alpha_num() - num_elec_diff
|
||||
nbeta = ezfio.get_electrons_elec_beta_num() - num_elec_diff
|
||||
|
||||
ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
|
||||
ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
|
||||
ezfio.set_electrons_elec_alpha_num(nalpha)
|
||||
ezfio.set_electrons_elec_beta_num( nbeta )
|
||||
|
||||
# Change all the array 'cause EZFIO
|
||||
# v_kl (v, l) => v_kl(l,v)
|
||||
|
|
@ -408,8 +429,8 @@ if __name__ == '__main__':
|
|||
|
||||
file_ = get_full_path(arguments['<file.out>'])
|
||||
|
||||
if arguments["--ezfio"]:
|
||||
ezfio_file = get_full_path(arguments["--ezfio"])
|
||||
if arguments["-o"]:
|
||||
ezfio_file = get_full_path(arguments["<ezfio_directory>"])
|
||||
else:
|
||||
ezfio_file = "{0}.ezfio".format(file_)
|
||||
|
||||
|
|
@ -421,3 +442,4 @@ if __name__ == '__main__':
|
|||
print file_, 'recognized as', str(res_file).split('.')[-1].split()[0]
|
||||
|
||||
write_ezfio(res_file, ezfio_file)
|
||||
os.system("qp_run save_ortho_mos "+ezfio_file)
|
||||
|
|
|
|||
|
|
@ -1,6 +1,10 @@
|
|||
open Qputils;;
|
||||
open Qptypes;;
|
||||
open Core.Std;;
|
||||
(*
|
||||
vim::syntax=ocaml
|
||||
*)
|
||||
|
||||
open Qputils
|
||||
open Qptypes
|
||||
open Core.Std
|
||||
|
||||
(** Interactive editing of the input.
|
||||
|
||||
|
|
@ -18,7 +22,7 @@ type keyword =
|
|||
| Mo_basis
|
||||
| Nuclei
|
||||
{keywords}
|
||||
;;
|
||||
|
||||
|
||||
|
||||
let keyword_to_string = function
|
||||
|
|
@ -28,7 +32,7 @@ let keyword_to_string = function
|
|||
| Mo_basis -> "MO basis"
|
||||
| Nuclei -> "Molecule"
|
||||
{keywords_to_string}
|
||||
;;
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -42,7 +46,7 @@ let file_header filename =
|
|||
Editing file `%s`
|
||||
|
||||
" filename
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Creates the header of a section *)
|
||||
|
|
@ -50,7 +54,7 @@ let make_header kw =
|
|||
let s = keyword_to_string kw in
|
||||
let l = String.length s in
|
||||
"\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Returns the rst string of section [s] *)
|
||||
|
|
@ -82,7 +86,7 @@ let get s =
|
|||
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
|
||||
in
|
||||
rst
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Applies the changes from the string [str] corresponding to section [s] *)
|
||||
|
|
@ -121,7 +125,7 @@ let set str s =
|
|||
| Ao_basis -> () (* TODO *)
|
||||
| Mo_basis -> () (* TODO *)
|
||||
end
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Creates the temporary file for interactive editing *)
|
||||
|
|
@ -135,11 +139,19 @@ let create_temp_file ezfio_filename fields =
|
|||
)
|
||||
end
|
||||
; temp_filename
|
||||
;;
|
||||
|
||||
|
||||
|
||||
|
||||
let run check_only ezfio_filename =
|
||||
|
||||
let run check_only ?ndet ?state ezfio_filename =
|
||||
|
||||
(* Set check_only if the arguments are not empty *)
|
||||
let check_only =
|
||||
match ndet, state with
|
||||
| None, None -> check_only
|
||||
| _ -> true
|
||||
in
|
||||
|
||||
(* Open EZFIO *)
|
||||
if (not (Sys.file_exists_exn ezfio_filename)) then
|
||||
|
|
@ -147,6 +159,19 @@ let run check_only ezfio_filename =
|
|||
|
||||
Ezfio.set_file ezfio_filename;
|
||||
|
||||
begin
|
||||
match ndet with
|
||||
| None -> ()
|
||||
| Some n -> Input.Determinants_by_hand.update_ndet (Det_number.of_int n)
|
||||
end;
|
||||
|
||||
begin
|
||||
match state with
|
||||
| None -> ()
|
||||
| Some n -> Input.Determinants_by_hand.extract_state (States_number.of_int n)
|
||||
end;
|
||||
|
||||
|
||||
(*
|
||||
let output = (file_header ezfio_filename) :: (
|
||||
List.map ~f:get [
|
||||
|
|
@ -196,7 +221,7 @@ let run check_only ezfio_filename =
|
|||
|
||||
(* Remove temp_file *)
|
||||
Sys.remove temp_filename
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Create a backup file in case of an exception *)
|
||||
|
|
@ -207,7 +232,7 @@ let create_backup ezfio_filename =
|
|||
"
|
||||
ezfio_filename ezfio_filename ezfio_filename
|
||||
|> Sys.command_exn
|
||||
;;
|
||||
|
||||
|
||||
|
||||
(** Restore the backup file when an exception occuprs *)
|
||||
|
|
@ -215,7 +240,7 @@ let restore_backup ezfio_filename =
|
|||
Printf.sprintf "tar -zxf %s/backup.tgz"
|
||||
ezfio_filename
|
||||
|> Sys.command_exn
|
||||
;;
|
||||
|
||||
|
||||
|
||||
let spec =
|
||||
|
|
@ -223,12 +248,12 @@ let spec =
|
|||
empty
|
||||
+> flag "-c" no_arg
|
||||
~doc:"Checks the input data"
|
||||
(*
|
||||
+> flag "o" (optional string)
|
||||
~doc:"Prints output data"
|
||||
*)
|
||||
+> flag "ndet" (optional int)
|
||||
~doc:"int Truncate the wavefunction to the target number of determinants"
|
||||
+> flag "state" (optional int)
|
||||
~doc:"int Pick the state as a new wavefunction."
|
||||
+> anon ("ezfio_file" %: string)
|
||||
;;
|
||||
|
||||
|
||||
let command =
|
||||
Command.basic
|
||||
|
|
@ -245,9 +270,9 @@ Edit input data
|
|||
with
|
||||
| _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
|
||||
*)
|
||||
(fun c ezfio_file () ->
|
||||
(fun c ndet state ezfio_file () ->
|
||||
try
|
||||
run c ezfio_file ;
|
||||
run c ?ndet ?state ezfio_file ;
|
||||
(* create_backup ezfio_file; *)
|
||||
with
|
||||
| Failure exc
|
||||
|
|
@ -268,12 +293,12 @@ Edit input data
|
|||
raise e
|
||||
end
|
||||
)
|
||||
;;
|
||||
|
||||
|
||||
let () =
|
||||
Command.run command;
|
||||
exit 0
|
||||
;;
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -2,11 +2,16 @@ use bitmasks
|
|||
|
||||
BEGIN_PROVIDER [ integer, N_int ]
|
||||
implicit none
|
||||
include 'Utils/constants.include.F'
|
||||
BEGIN_DOC
|
||||
! Number of 64-bit integers needed to represent determinants as binary strings
|
||||
END_DOC
|
||||
N_int = (mo_tot_num-1)/bit_kind_size + 1
|
||||
call write_int(6,N_int, 'N_int')
|
||||
call write_int(6,N_int, 'N_int')
|
||||
if (N_int > N_int_max) then
|
||||
stop 'N_int > N_int_max'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -355,7 +355,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
|||
write(iunit,'(A)') trim(write_buffer)
|
||||
write_buffer = ' Iter'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' Energy Residual'
|
||||
write_buffer = trim(write_buffer)//' Energy Residual'
|
||||
enddo
|
||||
write(iunit,'(A)') trim(write_buffer)
|
||||
write_buffer = '===== '
|
||||
|
|
@ -502,7 +502,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
|||
endif
|
||||
enddo
|
||||
|
||||
write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st)
|
||||
write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,E16.6))') iter, to_print(:,1:N_st)
|
||||
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
|
||||
if (converged) then
|
||||
exit
|
||||
|
|
|
|||
|
|
@ -413,7 +413,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
|
|||
endif
|
||||
enddo
|
||||
|
||||
write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(1:3,1:N_st)
|
||||
write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st)
|
||||
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
|
||||
do k=1,N_st
|
||||
if (residual_norm(k) > 1.e8) then
|
||||
|
|
@ -838,7 +838,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
|
|||
endif
|
||||
enddo
|
||||
|
||||
write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(1:3,1:N_st)
|
||||
write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st)
|
||||
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
|
||||
do k=1,N_st
|
||||
if (residual_norm(k) > 1.e8) then
|
||||
|
|
|
|||
|
|
@ -90,7 +90,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
|
|||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,1)
|
||||
do sh=1,shortcut(0,2)
|
||||
do i=shortcut(sh,2),shortcut(sh+1,2)-1
|
||||
org_i = sort_idx(i,2)
|
||||
|
|
@ -123,7 +123,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
|
|||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,1)
|
||||
do sh=1,shortcut(0,1)
|
||||
do sh2=1,shortcut(0,1)
|
||||
if (sh==sh2) cycle
|
||||
|
|
@ -342,7 +342,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_
|
|||
! istep = 1+ int(workload*target_workload_inv)
|
||||
istep = 1
|
||||
do blockb2=0, istep-1
|
||||
write(tmp_task,'(3(I9,X),''|'',X)') sh, blockb2, istep
|
||||
write(tmp_task,'(3(I9,1X),''|'',1X)') sh, blockb2, istep
|
||||
task = task//tmp_task
|
||||
ipos += 32
|
||||
if (ipos+32 > iposmax) then
|
||||
|
|
@ -444,7 +444,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
|||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,4)
|
||||
do sh=1,shortcut(0,2)
|
||||
do i=shortcut(sh,2),shortcut(sh+1,2)-1
|
||||
org_i = sort_idx(i,2)
|
||||
|
|
@ -477,7 +477,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
|||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,4)
|
||||
do sh=1,shortcut(0,1)
|
||||
do sh2=1,shortcut(0,1)
|
||||
if (sh==sh2) cycle
|
||||
|
|
@ -609,3 +609,352 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
|||
deallocate (shortcut, sort_idx, sorted, version, ut)
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
|
||||
!
|
||||
! n : number of determinants
|
||||
!
|
||||
! H_jj : array of <j|H|j>
|
||||
!
|
||||
! S2_jj : array of <j|S^2|j>
|
||||
END_DOC
|
||||
integer, intent(in) :: N_st,sze_8
|
||||
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
|
||||
|
||||
|
||||
PROVIDE ref_bitmask_energy
|
||||
|
||||
double precision :: hij, s2
|
||||
integer :: i,j
|
||||
integer :: k_a, k_b, l_a, l_b, m_a, m_b
|
||||
integer :: degree, istate
|
||||
integer :: krow, kcol, krow_b, kcol_b
|
||||
integer :: lrow, lcol
|
||||
integer :: mrow, mcol
|
||||
integer(bit_kind) :: spindet(N_int)
|
||||
integer(bit_kind) :: tmp_det(N_int,2)
|
||||
integer(bit_kind) :: tmp_det2(N_int,2)
|
||||
integer(bit_kind) :: tmp_det3(N_int,2)
|
||||
integer(bit_kind), allocatable :: buffer(:,:)
|
||||
double precision :: ck(N_st), cl(N_st), cm(N_st)
|
||||
integer :: n_singles, n_doubles
|
||||
integer, allocatable :: singles(:), doubles(:)
|
||||
integer, allocatable :: idx(:), idx0(:)
|
||||
logical, allocatable :: is_single_a(:)
|
||||
|
||||
allocate( buffer(N_int,N_det_alpha_unique), &
|
||||
singles(N_det_alpha_unique), doubles(N_det_alpha_unique), &
|
||||
is_single_a(N_det_alpha_unique), &
|
||||
idx(N_det_alpha_unique), idx0(N_det_alpha_unique) )
|
||||
|
||||
v_0 = 0.d0
|
||||
|
||||
do k_a=1,N_det-1
|
||||
|
||||
! Initial determinant is at k_a in alpha-major representation
|
||||
! -----------------------------------------------------------------------
|
||||
|
||||
krow = psi_bilinear_matrix_rows(k_a)
|
||||
kcol = psi_bilinear_matrix_columns(k_a)
|
||||
|
||||
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
|
||||
tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
|
||||
|
||||
! Initial determinant is at k_b in beta-major representation
|
||||
! ----------------------------------------------------------------------
|
||||
|
||||
k_b = psi_bilinear_matrix_order_reverse(k_a)
|
||||
|
||||
! Diagonal contribution
|
||||
! ---------------------
|
||||
|
||||
double precision, external :: diag_H_mat_elem
|
||||
|
||||
v_0(k_a,1:N_st) = v_0(k_a,1:N_st) + diag_H_mat_elem(tmp_det,N_int) * &
|
||||
psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
|
||||
|
||||
! Get all single and double alpha excitations
|
||||
! ===========================================
|
||||
|
||||
spindet(1:N_int) = tmp_det(1:N_int,1)
|
||||
|
||||
! Loop inside the beta column to gather all the connected alphas
|
||||
i=1
|
||||
l_a = k_a+1
|
||||
lcol = psi_bilinear_matrix_columns(l_a)
|
||||
do while ( (lcol == kcol).and.(l_a <= N_det) )
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow)
|
||||
idx(i) = lrow
|
||||
i=i+1
|
||||
l_a = l_a + 1
|
||||
lcol = psi_bilinear_matrix_columns(l_a)
|
||||
enddo
|
||||
i = i-1
|
||||
|
||||
call get_all_spin_singles_and_doubles( &
|
||||
buffer, idx, spindet, N_int, i, &
|
||||
singles, doubles, n_singles, n_doubles )
|
||||
|
||||
! Compute Hij for all alpha singles
|
||||
! ----------------------------------
|
||||
|
||||
l_a = k_a
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
tmp_det2(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
|
||||
do i=1,n_singles
|
||||
do while ( lrow < singles(i) )
|
||||
l_a = l_a+1
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
enddo
|
||||
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
|
||||
call i_H_j_mono_spin( tmp_det, tmp_det2, N_int, 1, hij)
|
||||
v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
|
||||
enddo
|
||||
|
||||
! Compute Hij for all alpha doubles
|
||||
! ----------------------------------
|
||||
|
||||
l_a = k_a
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
do i=1,n_doubles
|
||||
do while (lrow < doubles(i))
|
||||
l_a = l_a+1
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
enddo
|
||||
call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, doubles(i)), N_int, hij)
|
||||
v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
|
||||
enddo
|
||||
|
||||
|
||||
|
||||
! Get all single and double beta excitations
|
||||
! ===========================================
|
||||
|
||||
spindet(1:N_int) = tmp_det(1:N_int,2)
|
||||
|
||||
! Loop inside the alpha row to gather all the connected betas
|
||||
i=1
|
||||
l_b = k_b+1
|
||||
lrow = psi_bilinear_matrix_transp_rows(l_b)
|
||||
do while ( (lrow == krow).and.(l_b <= N_det) )
|
||||
lcol = psi_bilinear_matrix_transp_columns(l_b)
|
||||
buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol)
|
||||
idx(i) = lcol
|
||||
i=i+1
|
||||
l_b = l_b + 1
|
||||
lrow = psi_bilinear_matrix_transp_rows(l_b)
|
||||
enddo
|
||||
i = i-1
|
||||
|
||||
call get_all_spin_singles_and_doubles( &
|
||||
buffer, idx, spindet, N_int, i, &
|
||||
singles, doubles, n_singles, n_doubles )
|
||||
|
||||
! Compute Hij for all beta singles
|
||||
! ----------------------------------
|
||||
|
||||
l_b = k_b
|
||||
lcol = psi_bilinear_matrix_transp_columns(l_b)
|
||||
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
|
||||
do i=1,n_singles
|
||||
do while ( lcol < singles(i) )
|
||||
l_b = l_b+1
|
||||
lcol = psi_bilinear_matrix_transp_columns(l_b)
|
||||
enddo
|
||||
tmp_det2(1:N_int,2) = psi_det_beta_unique (1:N_int, lcol)
|
||||
l_a = psi_bilinear_matrix_transp_order(l_b)
|
||||
call i_H_j_mono_spin( tmp_det, tmp_det2, N_int, 2, hij)
|
||||
v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
|
||||
enddo
|
||||
|
||||
! Compute Hij for all beta doubles
|
||||
! ----------------------------------
|
||||
|
||||
l_b = k_b
|
||||
lcol = psi_bilinear_matrix_transp_columns(l_b)
|
||||
do i=1,n_doubles
|
||||
do while (lcol < doubles(i))
|
||||
l_b = l_b+1
|
||||
lcol = psi_bilinear_matrix_transp_columns(l_b)
|
||||
enddo
|
||||
l_a = psi_bilinear_matrix_transp_order(l_b)
|
||||
call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, doubles(i)), N_int, hij)
|
||||
v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
|
||||
enddo
|
||||
|
||||
end do
|
||||
|
||||
|
||||
! Alpha/Beta double excitations
|
||||
! =============================
|
||||
|
||||
do i=1,N_det_beta_unique
|
||||
idx0(i) = i
|
||||
enddo
|
||||
is_single_a(:) = .False.
|
||||
|
||||
k_a=1
|
||||
do i=1,N_det_beta_unique
|
||||
|
||||
! Select a beta determinant
|
||||
! -------------------------
|
||||
|
||||
spindet(1:N_int) = psi_det_beta_unique(1:N_int, i)
|
||||
tmp_det(1:N_int,2) = spindet(1:N_int)
|
||||
|
||||
call get_all_spin_singles( &
|
||||
psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, &
|
||||
singles, n_singles )
|
||||
|
||||
do j=1,n_singles
|
||||
is_single_a( singles(j) ) = .True.
|
||||
enddo
|
||||
|
||||
! For all alpha.beta pairs with the selected beta
|
||||
! -----------------------------------------------
|
||||
|
||||
kcol = psi_bilinear_matrix_columns(k_a)
|
||||
do while (kcol < i)
|
||||
k_a = k_a+1
|
||||
if (k_a > N_det) exit
|
||||
kcol = psi_bilinear_matrix_columns(k_a)
|
||||
enddo
|
||||
|
||||
do while (kcol == i)
|
||||
|
||||
krow = psi_bilinear_matrix_rows(k_a)
|
||||
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
|
||||
|
||||
! Loop over all alpha.beta pairs with a single exc alpha
|
||||
! ------------------------------------------------------
|
||||
|
||||
l_a = k_a+1
|
||||
if (l_a > N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
lcol = psi_bilinear_matrix_columns(l_a)
|
||||
|
||||
do while (lrow == krow)
|
||||
|
||||
! Loop over all alpha.beta pairs with a single exc alpha
|
||||
! ------------------------------------------------------
|
||||
if (is_single_a(lrow)) then
|
||||
|
||||
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow)
|
||||
|
||||
! Build list of singly excited beta
|
||||
! ---------------------------------
|
||||
|
||||
m_b = psi_bilinear_matrix_order_reverse(l_a)
|
||||
m_b = m_b+1
|
||||
j=1
|
||||
do while ( (mrow == lrow) )
|
||||
mcol = psi_bilinear_matrix_transp_columns(m_b)
|
||||
buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol)
|
||||
idx(j) = mcol
|
||||
j = j+1
|
||||
m_b = m_b+1
|
||||
if (m_b <= N_det) exit
|
||||
mrow = psi_bilinear_matrix_transp_rows(m_b)
|
||||
enddo
|
||||
j=j-1
|
||||
|
||||
call get_all_spin_singles( &
|
||||
buffer, idx, tmp_det(1,2), N_int, j, &
|
||||
doubles, n_doubles)
|
||||
|
||||
! Compute Hij for all doubles
|
||||
! ---------------------------
|
||||
|
||||
m_b = psi_bilinear_matrix_order(l_a)+1
|
||||
mcol = psi_bilinear_matrix_transp_columns(m_b)
|
||||
do j=1,n_doubles
|
||||
tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) )
|
||||
call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
|
||||
do while (mcol /= doubles(j))
|
||||
m_b = m_b+1
|
||||
if (m_b > N_det) exit
|
||||
mcol = psi_bilinear_matrix_transp_columns(m_b)
|
||||
enddo
|
||||
m_a = psi_bilinear_matrix_order_reverse(m_b)
|
||||
! v_0(m_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
|
||||
! v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(m_a,1:N_st)
|
||||
enddo
|
||||
|
||||
endif
|
||||
l_a = l_a+1
|
||||
if (l_a > N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l_a)
|
||||
lcol = psi_bilinear_matrix_columns(l_a)
|
||||
enddo
|
||||
|
||||
k_b = k_b+1
|
||||
if (k_b > N_det) exit
|
||||
kcol = psi_bilinear_matrix_transp_columns(k_b)
|
||||
enddo
|
||||
|
||||
do j=1,n_singles
|
||||
is_single_a( singles(j) ) = .False.
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: N_st,n,Nint, sze_8
|
||||
integer(bit_kind), intent(in) :: keys_tmp(Nint,2,n)
|
||||
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
|
||||
double precision, intent(in) :: u_0(sze_8,N_st)
|
||||
double precision, intent(in) :: H_jj(n), S2_jj(n)
|
||||
|
||||
PROVIDE ref_bitmask_energy
|
||||
|
||||
double precision, allocatable :: vt(:,:)
|
||||
integer, allocatable :: idx(:)
|
||||
integer :: i,j, jj
|
||||
double precision :: hij
|
||||
|
||||
do i=1,n
|
||||
v_0(i,:) = H_jj(i) * u_0(i,:)
|
||||
enddo
|
||||
|
||||
allocate(idx(0:n), vt(N_st,n))
|
||||
Vt = 0.d0
|
||||
do i=2,n
|
||||
idx(0) = i
|
||||
call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
|
||||
do jj=1,idx(0)
|
||||
j = idx(jj)
|
||||
double precision :: phase
|
||||
integer :: degree
|
||||
integer :: exc(0:2,2,2)
|
||||
call get_excitation(keys_tmp(1,1,j),keys_tmp(1,1,i),exc,degree,phase,Nint)
|
||||
if ((degree == 2).and.(exc(0,1,1)==1)) cycle
|
||||
! if (exc(0,1,2) /= 0) cycle
|
||||
call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
|
||||
vt (:,i) = vt (:,i) + hij*u_0(j,:)
|
||||
vt (:,j) = vt (:,j) + hij*u_0(i,:)
|
||||
enddo
|
||||
enddo
|
||||
do i=1,n
|
||||
v_0(i,:) = v_0(i,:) + vt(:,i)
|
||||
enddo
|
||||
end
|
||||
|
||||
|
|
|
|||
|
|
@ -78,96 +78,68 @@ END_PROVIDER
|
|||
double precision :: ck, cl, ckl
|
||||
double precision :: phase
|
||||
integer :: h1,h2,p1,p2,s1,s2, degree
|
||||
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int,2)
|
||||
integer :: exc(0:2,2,2),n_occ(2)
|
||||
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int)
|
||||
integer :: exc(0:2,2),n_occ(2)
|
||||
double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:)
|
||||
integer :: krow, kcol, lrow, lcol
|
||||
|
||||
one_body_dm_mo_alpha = 0.d0
|
||||
one_body_dm_mo_beta = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
|
||||
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
|
||||
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do k=1,N_det
|
||||
krow = psi_bilinear_matrix_rows(k)
|
||||
kcol = psi_bilinear_matrix_columns(k)
|
||||
tmp_det(:,1) = psi_det(:,1, krow)
|
||||
tmp_det(:,2) = psi_det(:,2, kcol)
|
||||
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
PROVIDE psi_det
|
||||
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
do while ( lcol == kcol )
|
||||
tmp_det2(:,1) = psi_det(:,1, lrow)
|
||||
tmp_det2(:,2) = psi_det(:,2, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
enddo
|
||||
one_body_dm_mo_alpha = 0.d0
|
||||
one_body_dm_mo_beta = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
|
||||
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
|
||||
!$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique, &
|
||||
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do k=1,N_det
|
||||
krow = psi_bilinear_matrix_rows(k)
|
||||
kcol = psi_bilinear_matrix_columns(k)
|
||||
tmp_det(:,1) = psi_det_alpha_unique(:,krow)
|
||||
tmp_det(:,2) = psi_det_beta_unique (:,kcol)
|
||||
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
do while ( lrow == krow )
|
||||
tmp_det2(:,1) = psi_det(:,1, lrow)
|
||||
tmp_det2(:,2) = psi_det(:,2, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
enddo
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
! Fix beta determinant, loop over alphas
|
||||
do while ( lcol == kcol )
|
||||
tmp_det2(:) = psi_det_alpha_unique(:, lrow)
|
||||
call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
exc = 0
|
||||
call get_mono_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
|
||||
call decode_exc_spin(exc,h1,p1,h2,p2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
|
|
@ -364,3 +336,74 @@ END_PROVIDER
|
|||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_old, (mo_tot_num_align,mo_tot_num,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_old, (mo_tot_num_align,mo_tot_num,N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha and beta one-body density matrix for each state
|
||||
END_DOC
|
||||
|
||||
integer :: j,k,l,m
|
||||
integer :: occ(N_int*bit_kind_size,2)
|
||||
double precision :: ck, cl, ckl
|
||||
double precision :: phase
|
||||
integer :: h1,h2,p1,p2,s1,s2, degree
|
||||
integer :: exc(0:2,2,2),n_occ(2)
|
||||
double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:)
|
||||
|
||||
one_body_dm_mo_alpha_old = 0.d0
|
||||
one_body_dm_mo_beta_old = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha_old,one_body_dm_mo_beta_old,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do k=1,N_det
|
||||
call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_coef(k,m)*psi_coef(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
do l=1,k-1
|
||||
call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
|
||||
if (degree /= 1) then
|
||||
cycle
|
||||
endif
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_coef(k,m) * psi_coef(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_alpha_old(:,:,:) = one_body_dm_mo_alpha_old(:,:,:) + tmp_a(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_beta_old(:,:,:) = one_body_dm_mo_beta_old(:,:,:) + tmp_b(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
deallocate(tmp_a,tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
END_PROVIDER
|
||||
|
|
|
|||
|
|
@ -23,7 +23,7 @@ BEGIN_PROVIDER [ integer, N_det ]
|
|||
! Number of determinants in the wave function
|
||||
END_DOC
|
||||
logical :: exists
|
||||
character*64 :: label
|
||||
character*(64) :: label
|
||||
PROVIDE ezfio_filename
|
||||
PROVIDE nproc
|
||||
if (read_wf) then
|
||||
|
|
@ -88,7 +88,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
|
|||
END_DOC
|
||||
integer :: i
|
||||
logical :: exists
|
||||
character*64 :: label
|
||||
character*(64) :: label
|
||||
|
||||
psi_det = 0_bit_kind
|
||||
if (read_wf) then
|
||||
|
|
|
|||
|
|
@ -252,8 +252,8 @@ 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
|
||||
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys)
|
||||
double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates)
|
||||
double precision, intent(out) :: s2(nstates,nstates)
|
||||
double precision :: s2_tmp,accu
|
||||
|
|
@ -344,7 +344,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nsta
|
|||
|
||||
print*,'S^2 matrix in the basis of the states considered'
|
||||
do i = 1, nstates
|
||||
write(*,'(100(F5.2,X))')s2(i,:)
|
||||
write(*,'(100(F5.2,1X))')s2(i,:)
|
||||
enddo
|
||||
|
||||
double precision :: accu_precision_diag,accu_precision_of_diag
|
||||
|
|
@ -370,7 +370,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nsta
|
|||
|
||||
print*,'Modified S^2 matrix that will be diagonalized'
|
||||
do i = 1, nstates
|
||||
write(*,'(10(F5.2,X))')s2(i,:)
|
||||
write(*,'(10(F5.2,1X))')s2(i,:)
|
||||
s2(i,i) = s2(i,i)
|
||||
enddo
|
||||
|
||||
|
|
|
|||
|
|
@ -1,32 +1,59 @@
|
|||
subroutine get_excitation_degree(key1,key2,degree,Nint)
|
||||
use bitmasks
|
||||
include 'Utils/constants.include.F'
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation degree between two determinants
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key2(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key1(Nint*2)
|
||||
integer(bit_kind), intent(in) :: key2(Nint*2)
|
||||
integer, intent(out) :: degree
|
||||
|
||||
integer(bit_kind) :: xorvec(2*N_int_max)
|
||||
integer :: l
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
degree = popcnt(xor( key1(1,1), key2(1,1))) + &
|
||||
popcnt(xor( key1(1,2), key2(1,2)))
|
||||
!DIR$ NOUNROLL
|
||||
do l=2,Nint
|
||||
degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + &
|
||||
popcnt(xor( key1(l,2), key2(l,2)))
|
||||
enddo
|
||||
ASSERT (degree >= 0)
|
||||
select case (Nint)
|
||||
|
||||
case (1)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
xorvec(2) = xor( key1(2), key2(2))
|
||||
degree = sum(popcnt(xorvec(1:2)))
|
||||
|
||||
case (2)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
xorvec(2) = xor( key1(2), key2(2))
|
||||
xorvec(3) = xor( key1(3), key2(3))
|
||||
xorvec(4) = xor( key1(4), key2(4))
|
||||
degree = sum(popcnt(xorvec(1:4)))
|
||||
|
||||
case (3)
|
||||
do l=1,6
|
||||
xorvec(l) = xor( key1(l), key2(l))
|
||||
enddo
|
||||
degree = sum(popcnt(xorvec(1:6)))
|
||||
|
||||
case (4)
|
||||
do l=1,8
|
||||
xorvec(l) = xor( key1(l), key2(l))
|
||||
enddo
|
||||
degree = sum(popcnt(xorvec(1:8)))
|
||||
|
||||
case default
|
||||
do l=1,ishft(Nint,1)
|
||||
xorvec(l) = xor( key1(l), key2(l))
|
||||
enddo
|
||||
degree = sum(popcnt(xorvec(1:l)))
|
||||
|
||||
end select
|
||||
|
||||
degree = ishft(degree,-1)
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
|
@ -139,72 +166,6 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|||
end
|
||||
|
||||
|
||||
subroutine decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Decodes the exc arrays returned by get_excitation.
|
||||
! h1,h2 : Holes
|
||||
! p1,p2 : Particles
|
||||
! s1,s2 : Spins (1:alpha, 2:beta)
|
||||
! degree : Degree of excitation
|
||||
END_DOC
|
||||
integer, intent(in) :: exc(0:2,2,2),degree
|
||||
integer*2, intent(out) :: h1,h2,p1,p2,s1,s2
|
||||
ASSERT (degree > 0)
|
||||
ASSERT (degree < 3)
|
||||
|
||||
select case(degree)
|
||||
case(2)
|
||||
if (exc(0,1,1) == 2) then
|
||||
h1 = exc(1,1,1)
|
||||
h2 = exc(2,1,1)
|
||||
p1 = exc(1,2,1)
|
||||
p2 = exc(2,2,1)
|
||||
s1 = 1
|
||||
s2 = 1
|
||||
else if (exc(0,1,2) == 2) then
|
||||
h1 = exc(1,1,2)
|
||||
h2 = exc(2,1,2)
|
||||
p1 = exc(1,2,2)
|
||||
p2 = exc(2,2,2)
|
||||
s1 = 2
|
||||
s2 = 2
|
||||
else
|
||||
h1 = exc(1,1,1)
|
||||
h2 = exc(1,1,2)
|
||||
p1 = exc(1,2,1)
|
||||
p2 = exc(1,2,2)
|
||||
s1 = 1
|
||||
s2 = 2
|
||||
endif
|
||||
case(1)
|
||||
if (exc(0,1,1) == 1) then
|
||||
h1 = exc(1,1,1)
|
||||
h2 = 0
|
||||
p1 = exc(1,2,1)
|
||||
p2 = 0
|
||||
s1 = 1
|
||||
s2 = 0
|
||||
else
|
||||
h1 = exc(1,1,2)
|
||||
h2 = 0
|
||||
p1 = exc(1,2,2)
|
||||
p2 = 0
|
||||
s1 = 2
|
||||
s2 = 0
|
||||
endif
|
||||
case(0)
|
||||
h1 = 0
|
||||
p1 = 0
|
||||
h2 = 0
|
||||
p2 = 0
|
||||
s1 = 0
|
||||
s2 = 0
|
||||
end select
|
||||
end
|
||||
|
||||
|
||||
subroutine get_double_excitation(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
|
@ -2154,8 +2115,8 @@ end
|
|||
subroutine get_phase(key1,key2,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2), key2(Nint,2)
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2), key2(Nint,2)
|
||||
double precision, intent(out) :: phase
|
||||
BEGIN_DOC
|
||||
! Returns the phase between key1 and key2
|
||||
|
|
@ -2224,3 +2185,423 @@ subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze)
|
|||
call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
|
||||
e_0 = u_dot_v(v_0,u_0,sze)
|
||||
end
|
||||
|
||||
|
||||
|
||||
! Spin-determinant routines
|
||||
! -------------------------
|
||||
|
||||
subroutine get_excitation_degree_spin(key1,key2,degree,Nint)
|
||||
use bitmasks
|
||||
include 'Utils/constants.include.F'
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation degree between two determinants
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key1(Nint)
|
||||
integer(bit_kind), intent(in) :: key2(Nint)
|
||||
integer, intent(out) :: degree
|
||||
|
||||
integer(bit_kind) :: xorvec(N_int_max)
|
||||
integer :: l
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
select case (Nint)
|
||||
|
||||
case (1)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
degree = popcnt(xorvec(1))
|
||||
|
||||
case (2)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
xorvec(2) = xor( key1(2), key2(2))
|
||||
degree = popcnt(xorvec(1))+popcnt(xorvec(2))
|
||||
|
||||
case (3)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
xorvec(2) = xor( key1(2), key2(2))
|
||||
xorvec(3) = xor( key1(3), key2(3))
|
||||
degree = sum(popcnt(xorvec(1:3)))
|
||||
|
||||
case (4)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
xorvec(2) = xor( key1(2), key2(2))
|
||||
xorvec(3) = xor( key1(3), key2(3))
|
||||
xorvec(4) = xor( key1(4), key2(4))
|
||||
degree = sum(popcnt(xorvec(1:4)))
|
||||
|
||||
case default
|
||||
do l=1,N_int
|
||||
xorvec(l) = xor( key1(l), key2(l))
|
||||
enddo
|
||||
degree = sum(popcnt(xorvec(1:Nint)))
|
||||
|
||||
end select
|
||||
|
||||
degree = ishft(degree,-1)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_excitation_spin(det1,det2,exc,degree,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation operators between two determinants and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
integer, intent(out) :: degree
|
||||
double precision, intent(out) :: phase
|
||||
! exc(number,hole/particle)
|
||||
! ex :
|
||||
! exc(0,1) = number of holes
|
||||
! exc(0,2) = number of particles
|
||||
! exc(1,2) = first particle
|
||||
! exc(1,1) = first hole
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call get_excitation_degree_spin(det1,det2,degree,Nint)
|
||||
select case (degree)
|
||||
|
||||
case (3:)
|
||||
degree = -1
|
||||
return
|
||||
|
||||
case (2)
|
||||
call get_double_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
return
|
||||
|
||||
case (1)
|
||||
call get_mono_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
return
|
||||
|
||||
case(0)
|
||||
return
|
||||
|
||||
end select
|
||||
end
|
||||
|
||||
subroutine decode_exc_spin(exc,h1,p1,h2,p2)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Decodes the exc arrays returned by get_excitation.
|
||||
! h1,h2 : Holes
|
||||
! p1,p2 : Particles
|
||||
END_DOC
|
||||
integer, intent(in) :: exc(0:2,2)
|
||||
integer, intent(out) :: h1,h2,p1,p2
|
||||
|
||||
select case (exc(0,1))
|
||||
case(2)
|
||||
h1 = exc(1,1)
|
||||
h2 = exc(2,1)
|
||||
p1 = exc(1,2)
|
||||
p2 = exc(2,2)
|
||||
case(1)
|
||||
h1 = exc(1,1)
|
||||
h2 = 0
|
||||
p1 = exc(1,2)
|
||||
p2 = 0
|
||||
case default
|
||||
h1 = 0
|
||||
p1 = 0
|
||||
h2 = 0
|
||||
p2 = 0
|
||||
end select
|
||||
end
|
||||
|
||||
|
||||
subroutine get_double_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the two excitation operators between two doubly excited spin-determinants
|
||||
! and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: tz
|
||||
integer :: l, idx_hole, idx_particle, ishift
|
||||
integer :: nperm
|
||||
integer :: i,j,k,m,n
|
||||
integer :: high, low
|
||||
integer :: a,b,c,d
|
||||
integer(bit_kind) :: hole, particle, tmp
|
||||
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
nperm = 0
|
||||
exc(0,1) = 0
|
||||
exc(0,2) = 0
|
||||
|
||||
idx_particle = 0
|
||||
idx_hole = 0
|
||||
ishift = 1-bit_kind_size
|
||||
do l=1,Nint
|
||||
ishift = ishift + bit_kind_size
|
||||
if (det1(l) == det2(l)) then
|
||||
cycle
|
||||
endif
|
||||
tmp = xor( det1(l), det2(l) )
|
||||
particle = iand(tmp, det2(l))
|
||||
hole = iand(tmp, det1(l))
|
||||
do while (particle /= 0_bit_kind)
|
||||
tz = trailz(particle)
|
||||
idx_particle = idx_particle + 1
|
||||
exc(0,2) = exc(0,2) + 1
|
||||
exc(idx_particle,2) = tz+ishift
|
||||
particle = iand(particle,particle-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1),exc(0,2))==2) then ! exc(0,1)==2 or exc(0,2)==2
|
||||
exit
|
||||
endif
|
||||
do while (hole /= 0_bit_kind)
|
||||
tz = trailz(hole)
|
||||
idx_hole = idx_hole + 1
|
||||
exc(0,1) = exc(0,1) + 1
|
||||
exc(idx_hole,1) = tz+ishift
|
||||
hole = iand(hole,hole-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1),exc(0,2))==2) then ! exc(0,1)==2 or exc(0,2)==2
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
select case (exc(0,1))
|
||||
|
||||
case(1)
|
||||
low = min(exc(1,1), exc(1,2))
|
||||
high = max(exc(1,1), exc(1,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
|
||||
if (j==k) then
|
||||
nperm = nperm + popcnt(iand(det1(j), &
|
||||
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
|
||||
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k), &
|
||||
ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
|
||||
endif
|
||||
do i=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(i))
|
||||
end do
|
||||
endif
|
||||
|
||||
case (2)
|
||||
|
||||
do i=1,2
|
||||
low = min(exc(i,1), exc(i,2))
|
||||
high = max(exc(i,1), exc(i,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
|
||||
if (j==k) then
|
||||
nperm = nperm + popcnt(iand(det1(j), &
|
||||
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
|
||||
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k), &
|
||||
ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
|
||||
endif
|
||||
do l=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(l))
|
||||
end do
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
a = min(exc(1,1), exc(1,2))
|
||||
b = max(exc(1,1), exc(1,2))
|
||||
c = min(exc(2,1), exc(2,2))
|
||||
d = max(exc(2,1), exc(2,2))
|
||||
if (c>a .and. c<b .and. d>b) then
|
||||
nperm = nperm + 1
|
||||
endif
|
||||
end select
|
||||
|
||||
phase = phase_dble(iand(nperm,1))
|
||||
|
||||
end
|
||||
|
||||
subroutine get_mono_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation operator between two singly excited determinants and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: tz
|
||||
integer :: l, idx_hole, idx_particle, ishift
|
||||
integer :: nperm
|
||||
integer :: i,j,k,m,n
|
||||
integer :: high, low
|
||||
integer :: a,b,c,d
|
||||
integer(bit_kind) :: hole, particle, tmp
|
||||
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
nperm = 0
|
||||
exc(0,1) = 0
|
||||
exc(0,2) = 0
|
||||
|
||||
ishift = 1-bit_kind_size
|
||||
do l=1,Nint
|
||||
ishift = ishift + bit_kind_size
|
||||
if (det1(l) == det2(l)) then
|
||||
cycle
|
||||
endif
|
||||
tmp = xor( det1(l), det2(l) )
|
||||
particle = iand(tmp, det2(l))
|
||||
hole = iand(tmp, det1(l))
|
||||
if (particle /= 0_bit_kind) then
|
||||
tz = trailz(particle)
|
||||
exc(0,2) = 1
|
||||
exc(1,2) = tz+ishift
|
||||
endif
|
||||
if (hole /= 0_bit_kind) then
|
||||
tz = trailz(hole)
|
||||
exc(0,1) = 1
|
||||
exc(1,1) = tz+ishift
|
||||
endif
|
||||
|
||||
if ( iand(exc(0,1),exc(0,2)) /= 1) then ! exc(0,1)/=1 and exc(0,2) /= 1
|
||||
cycle
|
||||
endif
|
||||
|
||||
low = min(exc(1,1),exc(1,2))
|
||||
high = max(exc(1,1),exc(1,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
if (j==k) then
|
||||
nperm = popcnt(iand(det1(j), &
|
||||
iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind)))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k),ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j),ibclr(-1_bit_kind,n)+1_bit_kind))
|
||||
endif
|
||||
do i=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(i))
|
||||
end do
|
||||
endif
|
||||
phase = phase_dble(iand(nperm,1))
|
||||
return
|
||||
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine i_H_j_mono_spin(key_i,key_j,Nint,spin,hij)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns <i|H|j> where i and j are determinants differing by a single excitation
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, spin
|
||||
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
|
||||
double precision, intent(out) :: hij
|
||||
|
||||
integer :: exc(0:2,2)
|
||||
double precision :: phase
|
||||
|
||||
PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
|
||||
|
||||
call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint)
|
||||
call get_mono_excitation_from_fock(key_i,key_j,exc(1,2),exc(1,1),spin,phase,hij)
|
||||
end
|
||||
|
||||
subroutine i_H_j_double_spin(key_i,key_j,Nint,hij)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns <i|H|j> where i and j are determinants differing by a same-spin double excitation
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint)
|
||||
double precision, intent(out) :: hij
|
||||
|
||||
integer :: exc(0:2,2)
|
||||
double precision :: phase
|
||||
double precision, external :: get_mo_bielec_integral
|
||||
|
||||
PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
|
||||
|
||||
call get_double_excitation_spin(key_i,key_j,exc,phase,Nint)
|
||||
hij = phase*(get_mo_bielec_integral( &
|
||||
exc(1,1), &
|
||||
exc(2,1), &
|
||||
exc(1,2), &
|
||||
exc(2,2), mo_integrals_map) - &
|
||||
get_mo_bielec_integral( &
|
||||
exc(1,1), &
|
||||
exc(2,1), &
|
||||
exc(2,2), &
|
||||
exc(1,2), mo_integrals_map) )
|
||||
end
|
||||
|
||||
subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns <i|H|j> where i and j are determinants differing by an opposite-spin double excitation
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
|
||||
double precision, intent(out) :: hij
|
||||
|
||||
integer :: exc(0:2,2,2)
|
||||
double precision :: phase
|
||||
double precision, external :: get_mo_bielec_integral
|
||||
|
||||
PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
|
||||
|
||||
call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint)
|
||||
call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase,Nint)
|
||||
if (exc(1,1,1) == exc(1,2,2)) then
|
||||
hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
|
||||
else if (exc(1,2,1) == exc(1,1,2)) then
|
||||
hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
|
||||
else
|
||||
hij = phase*get_mo_bielec_integral( &
|
||||
exc(1,1,1), &
|
||||
exc(1,1,2), &
|
||||
exc(1,2,1), &
|
||||
exc(1,2,2) ,mo_integrals_map)
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -386,8 +386,9 @@ END_PROVIDER
|
|||
!==============================================================================!
|
||||
|
||||
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows , (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order , (N_det) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
|
@ -395,19 +396,20 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
|
|||
! D_a^t C D_b
|
||||
!
|
||||
! Rows are alpha determinants and columns are beta.
|
||||
!
|
||||
! Order refers to psi_det
|
||||
END_DOC
|
||||
integer :: i,j,k, l
|
||||
integer(bit_kind) :: tmp_det(N_int,2)
|
||||
integer :: idx
|
||||
integer, external :: get_index_in_psi_det_sorted_bit
|
||||
|
||||
|
||||
PROVIDE psi_coef_sorted_bit
|
||||
|
||||
integer, allocatable :: iorder(:), to_sort(:)
|
||||
integer, allocatable :: to_sort(:)
|
||||
integer, external :: get_index_in_psi_det_alpha_unique
|
||||
integer, external :: get_index_in_psi_det_beta_unique
|
||||
allocate(iorder(N_det), to_sort(N_det))
|
||||
allocate(to_sort(N_det))
|
||||
do k=1,N_det
|
||||
i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int)
|
||||
j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int)
|
||||
|
|
@ -418,36 +420,41 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
|
|||
psi_bilinear_matrix_rows(k) = i
|
||||
psi_bilinear_matrix_columns(k) = j
|
||||
to_sort(k) = N_det_alpha_unique * (j-1) + i
|
||||
iorder(k) = k
|
||||
psi_bilinear_matrix_order(k) = k
|
||||
enddo
|
||||
call isort(to_sort, iorder, N_det)
|
||||
call iset_order(psi_bilinear_matrix_rows,iorder,N_det)
|
||||
call iset_order(psi_bilinear_matrix_columns,iorder,N_det)
|
||||
call isort(to_sort, psi_bilinear_matrix_order, N_det)
|
||||
call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det)
|
||||
call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det)
|
||||
do l=1,N_states
|
||||
call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det)
|
||||
call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det)
|
||||
enddo
|
||||
deallocate(iorder,to_sort)
|
||||
deallocate(to_sort)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows , (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order , (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
|
||||
! D_a^t C D_b
|
||||
!
|
||||
! Rows are Beta determinants and columns are alpha
|
||||
! Rows are Alpha determinants and columns are beta, but the matrix is stored in row major
|
||||
! format
|
||||
!
|
||||
! Order refers to psi_bilinear_matrix
|
||||
END_DOC
|
||||
integer :: i,j,k,l
|
||||
|
||||
|
||||
PROVIDE psi_coef_sorted_bit
|
||||
|
||||
integer, allocatable :: iorder(:), to_sort(:)
|
||||
allocate(iorder(N_det), to_sort(N_det))
|
||||
integer, allocatable :: to_sort(:)
|
||||
allocate(to_sort(N_det))
|
||||
do l=1,N_states
|
||||
do k=1,N_det
|
||||
psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l)
|
||||
|
|
@ -459,15 +466,18 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
|
|||
i = psi_bilinear_matrix_transp_columns(k)
|
||||
j = psi_bilinear_matrix_transp_rows (k)
|
||||
to_sort(k) = N_det_beta_unique * (j-1) + i
|
||||
iorder(k) = k
|
||||
psi_bilinear_matrix_transp_order(k) = k
|
||||
enddo
|
||||
call isort(to_sort, iorder, N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_rows,iorder,N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_columns,iorder,N_det)
|
||||
call isort(to_sort, psi_bilinear_matrix_transp_order, N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det)
|
||||
do l=1,N_states
|
||||
call dset_order(psi_bilinear_matrix_transp_values(1,l),iorder,N_det)
|
||||
call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det)
|
||||
enddo
|
||||
deallocate(iorder,to_sort)
|
||||
do k=1,N_det
|
||||
psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_transp_order(k)) = k
|
||||
enddo
|
||||
deallocate(to_sort)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
|
@ -552,7 +562,7 @@ subroutine generate_all_alpha_beta_det_products
|
|||
! Create a wave function from all possible alpha x beta determinants
|
||||
END_DOC
|
||||
integer :: i,j,k,l
|
||||
integer :: idx, iproc
|
||||
integer :: iproc
|
||||
integer, external :: get_index_in_psi_det_sorted_bit
|
||||
integer(bit_kind), allocatable :: tmp_det(:,:,:)
|
||||
logical, external :: is_in_wavefunction
|
||||
|
|
@ -561,7 +571,7 @@ subroutine generate_all_alpha_beta_det_products
|
|||
!$OMP PARALLEL DEFAULT(NONE) SHARED(psi_coef_sorted_bit,N_det_beta_unique,&
|
||||
!$OMP N_det_alpha_unique, N_int, psi_det_alpha_unique, psi_det_beta_unique,&
|
||||
!$OMP N_det) &
|
||||
!$OMP PRIVATE(i,j,k,l,tmp_det,idx,iproc)
|
||||
!$OMP PRIVATE(i,j,k,l,tmp_det,iproc)
|
||||
!$ iproc = omp_get_thread_num()
|
||||
allocate (tmp_det(N_int,2,N_det_alpha_unique))
|
||||
!$OMP DO
|
||||
|
|
@ -586,3 +596,782 @@ subroutine generate_all_alpha_beta_det_products
|
|||
end
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single and double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
! /!\ : The buffer is transposed !
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(Nint,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(Nint)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
select case (Nint)
|
||||
case (1)
|
||||
call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
return
|
||||
! case (2)
|
||||
! call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
! return
|
||||
case (3)
|
||||
call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
return
|
||||
end select
|
||||
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
|
||||
|
||||
do k=1,Nint
|
||||
do i=1,size_buffer
|
||||
xorvec(i, k) = xor( spindet(k), buffer(k,i) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
do k=2,Nint
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,k) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_singles(buffer, idx, spindet, Nint, size_buffer, singles, n_singles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(Nint,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(Nint)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
select case (Nint)
|
||||
case (1)
|
||||
call get_all_spin_singles_1(buffer, idx, spindet(1), size_buffer, singles, n_singles)
|
||||
return
|
||||
case (2)
|
||||
call get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
|
||||
return
|
||||
case (3)
|
||||
call get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
|
||||
return
|
||||
end select
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
|
||||
|
||||
do k=1,Nint
|
||||
do i=1,size_buffer
|
||||
xorvec(i, k) = xor( spindet(k), buffer(k,i) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
do k=2,Nint
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 2).and.(xorvec(i,k) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_doubles(buffer, idx, spindet, Nint, size_buffer, doubles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(Nint,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(Nint)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
select case (Nint)
|
||||
case (1)
|
||||
call get_all_spin_doubles_1(buffer, idx, spindet(1), size_buffer, doubles, n_doubles)
|
||||
return
|
||||
case (2)
|
||||
call get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
|
||||
return
|
||||
case (3)
|
||||
call get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
|
||||
return
|
||||
end select
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
|
||||
|
||||
do k=1,Nint
|
||||
do i=1,size_buffer
|
||||
xorvec(i, k) = xor( spindet(k), buffer(k,i) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
do k=2,Nint
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,k) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
enddo
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single and double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
! /!\ : The buffer is transposed !
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer
|
||||
integer, intent(in) :: idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:)
|
||||
integer :: degree
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i) = xor( spindet, buffer(i) )
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
n_doubles = 1
|
||||
|
||||
do i=1,size_buffer
|
||||
degree = popcnt(xorvec(i))
|
||||
if ( degree == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
if ( degree == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
n_doubles = n_doubles-1
|
||||
|
||||
deallocate(xorvec)
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_singles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:)
|
||||
|
||||
allocate( xorvec(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i) = xor( spindet, buffer(i) )
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
do i=1,size_buffer
|
||||
if ( popcnt(xorvec(i)) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:)
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
allocate( xorvec(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i) = xor( spindet, buffer(i) )
|
||||
enddo
|
||||
|
||||
n_doubles = 1
|
||||
|
||||
do i=1,size_buffer
|
||||
if ( popcnt(xorvec(i)) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
enddo
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single and double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
! /!\ : The buffer is transposed !
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(2,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(2)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(2,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(2)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 2).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(2,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(2)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
enddo
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
subroutine get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single and double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
! /!\ : The buffer is transposed !
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(3,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(3)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,3) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,3))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the single excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(3,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(3)
|
||||
integer, intent(out) :: singles(size_buffer)
|
||||
integer, intent(out) :: n_singles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 2).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 2).and.(xorvec(i,3) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,3))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_singles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 2 ) then
|
||||
singles(n_singles) = idx(i)
|
||||
n_singles = n_singles+1
|
||||
endif
|
||||
enddo
|
||||
n_singles = n_singles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Returns the indices of all the double excitations in the list of
|
||||
! unique alpha determinants.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: size_buffer, idx(size_buffer)
|
||||
integer(bit_kind), intent(in) :: buffer(3,size_buffer)
|
||||
integer(bit_kind), intent(in) :: spindet(3)
|
||||
integer, intent(out) :: doubles(size_buffer)
|
||||
integer, intent(out) :: n_doubles
|
||||
|
||||
integer :: i,k
|
||||
integer(bit_kind), allocatable :: xorvec(:,:)
|
||||
integer, allocatable :: degree(:)
|
||||
integer :: size_buffer_align
|
||||
|
||||
integer, external :: align_double
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
|
||||
|
||||
size_buffer_align = align_double(size_buffer)
|
||||
allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
|
||||
|
||||
do i=1,size_buffer
|
||||
xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
|
||||
xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
|
||||
xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if (xorvec(i,1) /= 0_8) then
|
||||
degree(i) = popcnt(xorvec(i,1))
|
||||
else
|
||||
degree(i) = 0
|
||||
endif
|
||||
enddo
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,2))
|
||||
endif
|
||||
enddo
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,size_buffer
|
||||
if ( (degree(i) <= 4).and.(xorvec(i,3) /= 0_8) ) then
|
||||
degree(i) = degree(i) + popcnt(xorvec(i,3))
|
||||
endif
|
||||
enddo
|
||||
|
||||
n_doubles = 1
|
||||
do i=1,size_buffer
|
||||
if ( degree(i) == 4 ) then
|
||||
doubles(n_doubles) = idx(i)
|
||||
n_doubles = n_doubles+1
|
||||
endif
|
||||
enddo
|
||||
n_doubles = n_doubles-1
|
||||
deallocate(xorvec)
|
||||
|
||||
end
|
||||
|
||||
|
|
|
|||
|
|
@ -2,7 +2,8 @@
|
|||
integer function n_open_shell(det_in,nint)
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer(bit_kind), intent(in) :: det_in(nint,2),nint
|
||||
integer, intent(in) :: nint
|
||||
integer(bit_kind), intent(in) :: det_in(nint,2)
|
||||
integer :: i
|
||||
n_open_shell = 0
|
||||
do i=1,Nint
|
||||
|
|
@ -13,7 +14,8 @@ end
|
|||
integer function n_closed_shell(det_in,nint)
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer(bit_kind), intent(in) :: det_in(nint,2),nint
|
||||
integer, intent(in) :: nint
|
||||
integer(bit_kind), intent(in) :: det_in(nint,2)
|
||||
integer :: i
|
||||
n_closed_shell = 0
|
||||
do i=1,Nint
|
||||
|
|
@ -24,7 +26,8 @@ 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, intent(in) :: nint
|
||||
integer(bit_kind), intent(in) :: det_in(nint,2)
|
||||
integer(bit_kind) :: det_tmp(nint,2)
|
||||
integer :: i
|
||||
n_closed_shell_cas = 0
|
||||
|
|
@ -44,8 +44,8 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
|
|||
l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
|
||||
i3 = i1 - ishft(i2*i2-i2,-1)
|
||||
k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
|
||||
j(1) = i2 - ishft(l(1)*l(1)-l(1),-1)
|
||||
i(1) = i3 - ishft(k(1)*k(1)-k(1),-1)
|
||||
j(1) = int(i2 - ishft(l(1)*l(1)-l(1),-1),4)
|
||||
i(1) = int(i3 - ishft(k(1)*k(1)-k(1),-1),4)
|
||||
|
||||
!ijkl
|
||||
i(2) = i(1) !ilkj
|
||||
|
|
|
|||
|
|
@ -51,7 +51,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
|
|||
print*, 'Providing the nuclear electron pseudo integrals (local)'
|
||||
|
||||
call wall_time(wall_1)
|
||||
wall_0 = wall_1
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
thread_num = 0
|
||||
|
|
@ -66,6 +65,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
|
|||
!$OMP wall_1)
|
||||
|
||||
!$ thread_num = omp_get_thread_num()
|
||||
|
||||
wall_0 = wall_1
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
|
||||
do j = 1, ao_num
|
||||
|
|
@ -148,7 +149,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
|
|||
print*, 'Providing the nuclear electron pseudo integrals (non-local)'
|
||||
|
||||
call wall_time(wall_1)
|
||||
wall_0 = wall_1
|
||||
call cpu_time(cpu_1)
|
||||
thread_num = 0
|
||||
|
||||
|
|
@ -164,6 +164,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
|
|||
|
||||
!$ thread_num = omp_get_thread_num()
|
||||
|
||||
wall_0 = wall_1
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
!
|
||||
do j = 1, ao_num
|
||||
|
|
|
|||
|
|
@ -42,7 +42,7 @@
|
|||
9;;
|
||||
END_TEMPLATE
|
||||
case default
|
||||
stop 'Error in ao_cart_to_sphe'
|
||||
stop 'Error in ao_cart_to_sphe : angular momentum too high'
|
||||
end select
|
||||
enddo
|
||||
|
||||
|
|
|
|||
|
|
@ -88,7 +88,7 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
|
|||
enddo
|
||||
endif
|
||||
do i=1,m
|
||||
write (output_mo_basis,'(I8,X,F16.10)') i,eigvalues(i)
|
||||
write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i)
|
||||
enddo
|
||||
write (output_mo_basis,'(A)') '======== ================'
|
||||
write (output_mo_basis,'(A)') ''
|
||||
|
|
@ -135,7 +135,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
|
|||
write (output_mo_basis,'(A)') '======== ================'
|
||||
|
||||
do i=1,m
|
||||
write (output_mo_basis,'(I8,X,F16.10)') i,D(i)
|
||||
write (output_mo_basis,'(I8,1X,F16.10)') i,D(i)
|
||||
enddo
|
||||
write (output_mo_basis,'(A)') '======== ================'
|
||||
write (output_mo_basis,'(A)') ''
|
||||
|
|
@ -215,7 +215,7 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,
|
|||
write (output_mo_basis,'(A)') ''
|
||||
write (output_mo_basis,'(A)') '======== ================'
|
||||
do i = 1, m
|
||||
write (output_mo_basis,'(I8,X,F16.10)') i,eigvalues(i)
|
||||
write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i)
|
||||
enddo
|
||||
write (output_mo_basis,'(A)') '======== ================'
|
||||
write (output_mo_basis,'(A)') ''
|
||||
|
|
|
|||
|
|
@ -37,8 +37,8 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num_aligned,3) ]
|
|||
enddo
|
||||
deallocate(buffer)
|
||||
|
||||
character*(64), parameter :: f = '(A16, 4(X,F12.6))'
|
||||
character*(64), parameter :: ft= '(A16, 4(X,A12 ))'
|
||||
character*(64), parameter :: f = '(A16, 4(1X,F12.6))'
|
||||
character*(64), parameter :: ft= '(A16, 4(1X,A12 ))'
|
||||
double precision, parameter :: a0= 0.529177249d0
|
||||
call write_time(output_Nuclei)
|
||||
write(output_Nuclei,'(A)') ''
|
||||
|
|
|
|||
|
|
@ -30,7 +30,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
|
|||
lwork = -1
|
||||
call dgesvd('A','A', m, n, A_tmp, LDA, &
|
||||
D, U, LDU, Vt, LDVt, work, lwork, info)
|
||||
lwork = work(1)
|
||||
lwork = int(work(1))
|
||||
deallocate(work)
|
||||
|
||||
allocate(work(lwork))
|
||||
|
|
@ -153,7 +153,7 @@ subroutine ortho_qr(A,LDA,m,n)
|
|||
allocate (jpvt(n), tau(n), work(1))
|
||||
LWORK=-1
|
||||
call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
|
||||
LWORK=2*WORK(1)
|
||||
LWORK=2*int(WORK(1))
|
||||
deallocate(WORK)
|
||||
allocate(WORK(LWORK))
|
||||
call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
|
||||
|
|
@ -297,7 +297,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
|
|||
print *, info, ': SVD failed'
|
||||
stop
|
||||
endif
|
||||
lwork = work(1)
|
||||
lwork = int(work(1))
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
integer, parameter :: max_dim = 511
|
||||
integer, parameter :: SIMD_vector = 32
|
||||
integer, parameter :: N_int_max = 16
|
||||
|
||||
double precision, parameter :: pi = dacos(-1.d0)
|
||||
double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))
|
||||
|
|
|
|||
|
|
@ -76,8 +76,8 @@ subroutine map_load_from_disk(filename,map)
|
|||
double precision :: x
|
||||
type(c_ptr) :: c_pointer(3)
|
||||
integer :: fd(3)
|
||||
integer*8 :: i,k, l
|
||||
integer :: n_elements, j
|
||||
integer*8 :: i,k,l
|
||||
integer*4 :: j,n_elements
|
||||
|
||||
|
||||
|
||||
|
|
@ -105,14 +105,14 @@ subroutine map_load_from_disk(filename,map)
|
|||
map % map(i) % value => map % consolidated_value ( map % consolidated_idx (i+1) :)
|
||||
map % map(i) % key => map % consolidated_key ( map % consolidated_idx (i+1) :)
|
||||
map % map(i) % sorted = .True.
|
||||
n_elements = map % consolidated_idx (i+2) - k
|
||||
n_elements = int( map % consolidated_idx (i+2) - k, 4)
|
||||
k = map % consolidated_idx (i+2)
|
||||
map % map(i) % map_size = n_elements
|
||||
map % map(i) % n_elements = n_elements
|
||||
! Load memory from disk
|
||||
do j=1,n_elements
|
||||
x = x + map % map(i) % value(j)
|
||||
l = iand(l,map % map(i) % key(j))
|
||||
l = iand(l,int(map % map(i) % key(j),8))
|
||||
if (map % map(i) % value(j) > 1.e30) then
|
||||
stop 'Error in integrals file'
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -53,17 +53,17 @@ module map_module
|
|||
end module map_module
|
||||
|
||||
|
||||
real function map_mb(map)
|
||||
double precision function map_mb(map)
|
||||
use map_module
|
||||
use omp_lib
|
||||
implicit none
|
||||
type (map_type), intent(in) :: map
|
||||
integer(map_size_kind) :: i
|
||||
|
||||
map_mb = 8+map_size_kind+map_size_kind+omp_lock_kind+4
|
||||
map_mb = dble(8+map_size_kind+map_size_kind+omp_lock_kind+4)
|
||||
do i=0,map%map_size
|
||||
map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) +&
|
||||
8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind
|
||||
map_mb = map_mb + dble(map%map(i)%map_size*(cache_key_kind+integral_kind) +&
|
||||
8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind)
|
||||
enddo
|
||||
map_mb = map_mb / (1024.d0*1024.d0)
|
||||
end
|
||||
|
|
@ -406,8 +406,8 @@ subroutine map_update(map, key, value, sze, thr)
|
|||
call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
|
||||
call cache_map_shrink(local_map,thr)
|
||||
endif
|
||||
cache_key = iand(key(i),map_mask)
|
||||
local_map%n_elements = local_map%n_elements + 1_8
|
||||
cache_key = int(iand(key(i),map_mask),2)
|
||||
local_map%n_elements = local_map%n_elements + 1
|
||||
local_map%value(local_map%n_elements) = value(i)
|
||||
local_map%key(local_map%n_elements) = cache_key
|
||||
local_map%sorted = .False.
|
||||
|
|
@ -464,7 +464,7 @@ subroutine map_append(map, key, value, sze)
|
|||
if (n_elements == map%map(idx_cache)%map_size) then
|
||||
call cache_map_reallocate(map%map(idx_cache), n_elements+ ishft(n_elements,-1))
|
||||
endif
|
||||
cache_key = iand(key(i),map_mask)
|
||||
cache_key = int(iand(key(i),map_mask),2)
|
||||
map%map(idx_cache)%value(n_elements) = value(i)
|
||||
map%map(idx_cache)%key(n_elements) = cache_key
|
||||
map%map(idx_cache)%n_elements = n_elements
|
||||
|
|
@ -615,7 +615,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
|
|||
idx = -1
|
||||
return
|
||||
endif
|
||||
cache_key = iand(key,map_mask)
|
||||
cache_key = int(iand(key,map_mask),2)
|
||||
ibegin = min(ibegin_in,sze)
|
||||
iend = min(iend_in,sze)
|
||||
if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
|
||||
|
|
@ -723,7 +723,7 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
|
|||
value = 0.d0
|
||||
return
|
||||
endif
|
||||
cache_key = iand(key,map_mask)
|
||||
cache_key = int(iand(key,map_mask),2)
|
||||
ibegin = min(ibegin_in,sze)
|
||||
iend = min(iend_in,sze)
|
||||
if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
|
||||
|
|
|
|||
|
|
@ -292,18 +292,17 @@ BEGIN_TEMPLATE
|
|||
! contains the new order of the elements.
|
||||
! iradix should be -1 in input.
|
||||
END_DOC
|
||||
$int_type, intent(in) :: isize
|
||||
$int_type, intent(inout) :: iorder(isize)
|
||||
$type, intent(inout) :: x(isize)
|
||||
integer*$int_type, intent(in) :: isize
|
||||
integer*$int_type, intent(inout) :: iorder(isize)
|
||||
integer*$type, intent(inout) :: x(isize)
|
||||
integer, intent(in) :: iradix
|
||||
integer :: iradix_new
|
||||
$type, allocatable :: x2(:), x1(:)
|
||||
$type :: i4
|
||||
$int_type, allocatable :: iorder1(:),iorder2(:)
|
||||
$int_type :: i0, i1, i2, i3, i
|
||||
integer*$type, allocatable :: x2(:), x1(:)
|
||||
integer*$type :: i4
|
||||
integer*$int_type, allocatable :: iorder1(:),iorder2(:)
|
||||
integer*$int_type :: i0, i1, i2, i3, i
|
||||
integer, parameter :: integer_size=$octets
|
||||
$type, parameter :: zero=$zero
|
||||
$type :: mask
|
||||
integer*$type :: mask
|
||||
integer :: nthreads, omp_get_num_threads
|
||||
!DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1
|
||||
|
||||
|
|
@ -311,16 +310,16 @@ BEGIN_TEMPLATE
|
|||
|
||||
! Find most significant bit
|
||||
|
||||
i0 = 0_8
|
||||
i4 = -1_8
|
||||
i0 = 0_$int_type
|
||||
i4 = -1_$type
|
||||
|
||||
do i=1,isize
|
||||
i4 = max(i4,x(i))
|
||||
enddo
|
||||
i3 = i4 ! Type conversion
|
||||
i3 = int(i4,$int_type)
|
||||
|
||||
iradix_new = integer_size-1-leadz(i3)
|
||||
mask = ibset(zero,iradix_new)
|
||||
mask = ibset(0_$type,iradix_new)
|
||||
nthreads = 1
|
||||
! nthreads = 1+ishft(omp_get_num_threads(),-1)
|
||||
|
||||
|
|
@ -331,22 +330,22 @@ BEGIN_TEMPLATE
|
|||
stop
|
||||
endif
|
||||
|
||||
i1=1_8
|
||||
i2=1_8
|
||||
i1=1_$int_type
|
||||
i2=1_$int_type
|
||||
|
||||
do i=1,isize
|
||||
if (iand(mask,x(i)) == zero) then
|
||||
if (iand(mask,x(i)) == 0_$type) then
|
||||
iorder1(i1) = iorder(i)
|
||||
x1(i1) = x(i)
|
||||
i1 = i1+1_8
|
||||
i1 = i1+1_$int_type
|
||||
else
|
||||
iorder2(i2) = iorder(i)
|
||||
x2(i2) = x(i)
|
||||
i2 = i2+1_8
|
||||
i2 = i2+1_$int_type
|
||||
endif
|
||||
enddo
|
||||
i1=i1-1_8
|
||||
i2=i2-1_8
|
||||
i1=i1-1_$int_type
|
||||
i2=i2-1_$int_type
|
||||
|
||||
do i=1,i1
|
||||
iorder(i0+i) = iorder1(i)
|
||||
|
|
@ -399,12 +398,12 @@ BEGIN_TEMPLATE
|
|||
endif
|
||||
|
||||
|
||||
mask = ibset(zero,iradix)
|
||||
mask = ibset(0_$type,iradix)
|
||||
i0=1
|
||||
i1=1
|
||||
|
||||
do i=1,isize
|
||||
if (iand(mask,x(i)) == zero) then
|
||||
if (iand(mask,x(i)) == 0_$type) then
|
||||
iorder(i0) = iorder(i)
|
||||
x(i0) = x(i)
|
||||
i0 = i0+1
|
||||
|
|
@ -443,12 +442,12 @@ BEGIN_TEMPLATE
|
|||
|
||||
end
|
||||
|
||||
SUBST [ X, type, octets, is_big, big, int_type, zero ]
|
||||
i ; integer ; 32 ; .False. ; ; integer ; 0;;
|
||||
i8 ; integer*8 ; 32 ; .False. ; ; integer ; 0_8;;
|
||||
i2 ; integer*2 ; 32 ; .False. ; ; integer ; 0;;
|
||||
i ; integer ; 64 ; .True. ; _big ; integer*8 ; 0 ;;
|
||||
i8 ; integer*8 ; 64 ; .True. ; _big ; integer*8 ; 0_8 ;;
|
||||
SUBST [ X, type, octets, is_big, big, int_type ]
|
||||
i ; 4 ; 32 ; .False. ; ; 4 ;;
|
||||
i8 ; 8 ; 32 ; .False. ; ; 4 ;;
|
||||
i2 ; 2 ; 32 ; .False. ; ; 4 ;;
|
||||
i ; 4 ; 64 ; .True. ; _big ; 8 ;;
|
||||
i8 ; 8 ; 64 ; .True. ; _big ; 8 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,11 +1,8 @@
|
|||
use f77_zmq
|
||||
use omp_lib
|
||||
|
||||
integer, pointer :: thread_id
|
||||
integer(omp_lock_kind) :: zmq_lock
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
|
||||
&BEGIN_PROVIDER [ integer(omp_lock_kind), zmq_lock ]
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
|
@ -407,7 +404,9 @@ subroutine end_zmq_sub_socket(zmq_socket_sub)
|
|||
integer(ZMQ_PTR), intent(in) :: zmq_socket_sub
|
||||
integer :: rc
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
rc = f77_zmq_close(zmq_socket_sub)
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (rc /= 0) then
|
||||
print *, 'f77_zmq_close(zmq_socket_sub)'
|
||||
stop 'error'
|
||||
|
|
@ -426,7 +425,9 @@ subroutine end_zmq_pair_socket(zmq_socket_pair)
|
|||
integer :: rc
|
||||
character*(8), external :: zmq_port
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
rc = f77_zmq_close(zmq_socket_pair)
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (rc /= 0) then
|
||||
print *, 'f77_zmq_close(zmq_socket_pair)'
|
||||
stop 'error'
|
||||
|
|
@ -444,7 +445,9 @@ subroutine end_zmq_pull_socket(zmq_socket_pull)
|
|||
integer :: rc
|
||||
character*(8), external :: zmq_port
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
rc = f77_zmq_close(zmq_socket_pull)
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (rc /= 0) then
|
||||
print *, 'f77_zmq_close(zmq_socket_pull)'
|
||||
stop 'error'
|
||||
|
|
@ -469,7 +472,9 @@ subroutine end_zmq_push_socket(zmq_socket_push,thread)
|
|||
stop 'Unable to set ZMQ_LINGER on push socket'
|
||||
endif
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
rc = f77_zmq_close(zmq_socket_push)
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (rc /= 0) then
|
||||
print *, 'f77_zmq_close(zmq_socket_push)'
|
||||
stop 'error'
|
||||
|
|
@ -500,10 +505,17 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,name_in)
|
|||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
zmq_context = f77_zmq_ctx_new ()
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (zmq_context == 0_ZMQ_PTR) then
|
||||
stop 'ZMQ_PTR is null'
|
||||
endif
|
||||
! rc = f77_zmq_ctx_set(zmq_context, ZMQ_IO_THREADS, nproc)
|
||||
! if (rc /= 0) then
|
||||
! print *, 'Unable to set the number of ZMQ IO threads to', nproc
|
||||
! endif
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
name = name_in
|
||||
sze = len(trim(name))
|
||||
|
|
@ -584,7 +596,10 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,name_in)
|
|||
zmq_state = 'No_state'
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
|
||||
call omp_set_lock(zmq_lock)
|
||||
rc = f77_zmq_ctx_term(zmq_context)
|
||||
zmq_context = 0_ZMQ_PTR
|
||||
call omp_unset_lock(zmq_lock)
|
||||
if (rc /= 0) then
|
||||
print *, 'Unable to terminate ZMQ context'
|
||||
stop 'error'
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
|
|||
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]"
|
||||
qp_run cassd_zmq $INPUT
|
||||
energy="$(ezfio get cas_sd_zmq energy_pt2)"
|
||||
eq $energy -76.231084536315 5.E-5
|
||||
eq $energy -76.231248286858 5.E-5
|
||||
|
||||
ezfio set determinants n_det_max 1024
|
||||
ezfio set determinants read_wf True
|
||||
|
|
@ -21,6 +21,6 @@ source $QP_ROOT/tests/bats/common.bats.sh
|
|||
qp_run cassd_zmq $INPUT
|
||||
ezfio set determinants read_wf False
|
||||
energy="$(ezfio get cas_sd_zmq energy)"
|
||||
eq $energy -76.2225863580749 2.E-5
|
||||
eq $energy -76.2225678834779 2.E-5
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -42,11 +42,13 @@ function run_FCI_ZMQ() {
|
|||
qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
|
||||
}
|
||||
@test "FCI H2O cc-pVDZ" {
|
||||
run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
|
||||
run_FCI h2o.ezfio 2000 -76.1253758241716 -76.1258130146102
|
||||
}
|
||||
|
||||
|
||||
|
||||
@test "FCI-ZMQ H2O cc-pVDZ" {
|
||||
run_FCI_ZMQ h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
|
||||
run_FCI_ZMQ h2o.ezfio 2000 -76.1250552686394 -76.1258817228809
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
|
|||
ezfio set mrcepa0 n_it_max_dressed_ci 3
|
||||
qp_run $EXE $INPUT
|
||||
energy="$(ezfio get mrcepa0 energy_pt2)"
|
||||
eq $energy -76.2381673136696 2.e-4
|
||||
eq $energy -76.2381754078899 1.e-4
|
||||
}
|
||||
|
||||
@test "MRSC2 H2O cc-pVDZ" {
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user