LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Merge branch 'scemama-master'

Browse Source
This commit is contained in:
Emmanuel Giner 2017-04-13 16:27:35 +02:00
commit 72ed06bf19
57 changed files with 3459 additions and 542 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
config/gfortran_avx.cfg
View File

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

2
config/gfortran_debug.cfg
View File

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

86
ocaml/Input_determinants_by_hand.ml
View File

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

2
ocaml/Symmetry.ml
View File

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

8
ocaml/TaskServer.ml
View File

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

51
ocaml/qp_create_ezfio_from_xyz.ml
View File

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

8
plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
View File

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

36
plugins/CAS_SD_ZMQ/selection.irp.f
View File

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

28
plugins/DFT_Utils/grid_density.irp.f
View File

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

6
plugins/DFT_Utils/integration_radial.irp.f
View File

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

2
plugins/Full_CI_ZMQ/energy.irp.f
View File

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

27
plugins/Full_CI_ZMQ/pt2_stoch.irp.f
View File

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

44
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

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

18
plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
View File

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

1
plugins/Full_CI_ZMQ/run_selection_slave.irp.f
View File

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

1112
plugins/Full_CI_ZMQ/selection.irp.f
View File

File diff suppressed because it is too large Load Diff

22
plugins/Full_CI_ZMQ/selection_buffer.irp.f
View File

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

6
plugins/Full_CI_ZMQ/selection_types.f90
View File

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

27
plugins/Full_CI_ZMQ/zmq_selection.irp.f
View File

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

78
plugins/MRCC_Utils/amplitudes.irp.f
View File

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

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

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

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

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

2
plugins/MRPT_Utils/excitations_cas.irp.f
View File

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

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

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

2
plugins/Properties/delta_rho.irp.f
View File

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

2
plugins/Properties/hyperfine_constants.irp.f
View File

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

4
plugins/Properties/mulliken.irp.f
View File

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

6
plugins/analyze_wf/analyze_wf.irp.f
View File

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

2
plugins/mrcc_selected/ezfio_interface.irp.f
View File

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

2
plugins/mrcepa0/EZFIO.cfg
View File

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

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

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

46
scripts/ezfio_interface/qp_convert_output_to_ezfio.py
View File

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

71
scripts/ezfio_interface/qp_edit_template
View 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
;;

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

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

4
src/Davidson/diagonalization.irp.f
View File

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

4
src/Davidson/diagonalization_hs2.irp.f
View File

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

359
src/Davidson/u0Hu0.irp.f
View File

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

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

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

4
src/Determinants/determinants.irp.f
View File

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

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

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

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

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

833
src/Determinants/spindeterminants.irp.f
View File

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

9
src/Determinants/usefull_for_ovb.irp.f → src/Determinants/useful_for_ovb.irp.f
View File

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

4
src/Integrals_Bielec/map_integrals.irp.f
View File

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

5
src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
View File

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

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

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

6
src/MO_Basis/utils.irp.f
View File

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

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

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

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

@ -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
src/Utils/constants.include.F
View File

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

8
src/Utils/map_functions.irp.f
View File

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

18
src/Utils/map_module.f90
View File

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

55
src/Utils/sort.irp.f
View File

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

25
src/ZMQ/utils.irp.f
View File

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

4
tests/bats/cassd.bats
View File

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

6
tests/bats/fci.bats
View File

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

2
tests/bats/mrcepa0.bats
View File

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