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

med

Browse Source
This commit is contained in:
Emmanuel Giner 2017-04-14 19:10:18 +02:00
commit 63b59f7d30
75 changed files with 11148 additions and 790 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

4
plugins/DFT_Utils/EZFIO.cfg
View File

@ -1,4 +0,0 @@
[energy]
type: double precision
doc: Calculated energy
interface: ezfio

6951
plugins/DFT_Utils/angular.f Normal file
View File

File diff suppressed because it is too large Load Diff

25
plugins/DFT_Utils/functional.irp.f Normal file
View File

@ -0,0 +1,25 @@
double precision function ex_lda(rho)
include 'constants.include.F'
implicit none
double precision, intent(in) :: rho
ex_lda = cst_lda * rho**(c_4_3)
end
BEGIN_PROVIDER [double precision, lda_exchange, (N_states)]
implicit none
integer :: i,j,k,l
double precision :: ex_lda
do l = 1, N_states
lda_exchange(l) = 0.d0
do j = 1, nucl_num
do i = 1, n_points_radial_grid
do k = 1, n_points_integration_angular
lda_exchange(l) += final_weight_functions_at_grid_points(k,i,j) * &
(ex_lda(one_body_dm_mo_alpha_at_grid_points(k,i,j,l)) + ex_lda(one_body_dm_mo_beta_at_grid_points(k,i,j,l)))
enddo
enddo
enddo
enddo
END_PROVIDER

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

@ -1,7 +1,7 @@
BEGIN_PROVIDER [integer, n_points_angular_grid]
BEGIN_PROVIDER [integer, n_points_integration_angular]
implicit none
n_points_angular_grid = 50
END_PROVIDER
n_points_integration_angular = 110
END_PROVIDER
BEGIN_PROVIDER [integer, n_points_radial_grid]
implicit none
@ -9,34 +9,52 @@ 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,3) ]
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)]
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)
angular_quadrature_points = 0.d0
weights_angular_points = 0.d0
!call cal_quad(n_points_integration_angular, angular_quadrature_points,weights_angular_points)
include 'constants.include.F'
integer :: i
integer :: i,n
double precision :: accu
double precision :: degre_rad
!degre_rad = 180.d0/pi
!accu = 0.d0
!do i = 1, n_points_integration_angular_lebedev
degre_rad = pi/180.d0
accu = 0.d0
double precision :: x(n_points_integration_angular),y(n_points_integration_angular),z(n_points_integration_angular),w(n_points_integration_angular)
call LD0110(X,Y,Z,W,N)
do i = 1, n_points_integration_angular
angular_quadrature_points(i,1) = x(i)
angular_quadrature_points(i,2) = y(i)
angular_quadrature_points(i,3) = z(i)
weights_angular_points(i) = w(i) * 4.d0 * pi
accu += w(i)
enddo
!do i = 1, n_points_integration_angular
! accu += weights_angular_integration_lebedev(i)
! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 2.d0 * pi
! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 4.d0 * pi
! angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) &
! * dsin ( degre_rad * phi_angular_integration_lebedev(i))
! angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) &
! * dsin ( degre_rad * phi_angular_integration_lebedev(i))
! angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i))
!!weights_angular_points(i) = weights_angular_integration_lebedev(i)
!!angular_quadrature_points(i,1) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) &
!! * dsin ( degre_rad * theta_angular_integration_lebedev(i))
!!angular_quadrature_points(i,2) = dsin ( degre_rad * phi_angular_integration_lebedev(i)) &
!! * dsin ( degre_rad * theta_angular_integration_lebedev(i))
!!angular_quadrature_points(i,3) = dcos ( degre_rad * theta_angular_integration_lebedev(i))
!enddo
!print*,'ANGULAR'
!print*,''
!print*,'accu = ',accu
!ASSERT( dabs(accu - 1.D0) < 1.d-10)
print*,'ANGULAR'
print*,''
print*,'accu = ',accu
ASSERT( dabs(accu - 1.D0) < 1.d-10)
END_PROVIDER
@ -63,7 +81,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,n_points_radial_grid,nucl_num)]
BEGIN_DOC
! points for integration over space
END_DOC
@ -79,7 +97,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 ! 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 +106,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,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 +120,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 ! 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)
@ -115,7 +133,6 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
enddo
accu = 1.d0/accu
weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu
! print*,weight_functions_at_grid_points(l,k,j)
enddo
enddo
enddo
@ -123,43 +140,64 @@ 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, final_weight_functions_at_grid_points, (n_points_integration_angular,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
! and the points are labelled by the other dimensions
END_DOC
implicit none
integer :: i,j,k,l,m
double precision :: r(3)
double precision :: accu,cell_function_becke
double precision :: tmp_array(nucl_num)
double precision :: contrib_integration,x
double precision :: derivative_knowles_function,knowles_function
! run over all points in space
do j = 1, nucl_num ! that are referred to each atom
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x) &
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
final_weight_functions_at_grid_points(k,i,j) = weights_angular_points(k) * weight_functions_at_grid_points(k,i,j) * contrib_integration * dr_radial_integral
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
implicit none
integer :: i,j,k,l,m,i_state
double precision :: contrib
double precision :: r(3)
double precision :: aos_array(ao_num),mos_array(mo_tot_num)
do i_state = 1, N_states
do j = 1, nucl_num
do k = 1, n_points_radial_grid -1
do l = 1, n_points_angular_grid
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
do k = 1, n_points_radial_grid
do l = 1, n_points_integration_angular
one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) = 0.d0
one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) = 0.d0
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)
! call give_all_aos_at_r(r,aos_array)
! do i = 1, ao_num
! do m = 1, ao_num
! contrib = aos_array(i) * aos_array(m)
! one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_ao_alpha(i,m) * contrib
! one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_ao_beta(i,m) * contrib
! enddo
! enddo
call give_all_mos_at_r(r,mos_array)
do i = 1, mo_tot_num
do m = 1, mo_tot_num
do m = 1, mo_tot_num
do i = 1, mo_tot_num
contrib = mos_array(i) * mos_array(m)
one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_mo_alpha(i,m) * contrib
one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_mo_beta(i,m) * contrib
one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) += one_body_dm_mo_alpha(i,m,i_state) * contrib
one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) += one_body_dm_mo_beta(i,m,i_state) * contrib
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER

9
plugins/DFT_Utils/integration_3d.irp.f
View File

@ -4,18 +4,11 @@ double precision function step_function_becke(x)
double precision :: f_function_becke
integer :: i,n_max_becke
!if(x.lt.-1.d0)then
! step_function_becke = 0.d0
!else if (x .gt.1)then
! step_function_becke = 0.d0
!else
step_function_becke = f_function_becke(x)
!!n_max_becke = 1
do i = 1, 4
do i = 1,5
step_function_becke = f_function_becke(step_function_becke)
enddo
step_function_becke = 0.5d0*(1.d0 - step_function_becke)
!endif
end
double precision function f_function_becke(x)

11
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), weights(n_points_integration_angular)
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 (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,14 +20,13 @@
! 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
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)
do k = 1, n_points_integration_angular
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j,1) * 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,1) * weight_functions_at_grid_points(k,i,j)
enddo
!
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
double precision :: contrib_integration
! print*,m_knowles
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x) &
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
integral_density_alpha_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_alpha

54
plugins/DFT_Utils/test_integration_3d_density.irp.f
View File

@ -4,13 +4,55 @@ program pouet
touch read_wf
print*,'m_knowles = ',m_knowles
call routine
call routine3
end
subroutine routine3
implicit none
integer :: i,j,k,l
double precision :: accu
accu = 0.d0
do j = 1, nucl_num ! that are referred to each atom
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
accu += final_weight_functions_at_grid_points(k,i,j) * one_body_dm_mo_alpha_at_grid_points(k,i,j,1)
enddo
enddo
enddo
print*, accu
print*, 'lda_exchange',lda_exchange
end
subroutine routine2
implicit none
integer :: i,j,k,l
double precision :: x,y,z
double precision :: r
double precision :: accu
accu = 0.d0
r = 1.d0
do k = 1, n_points_integration_angular
x = angular_quadrature_points(k,1) * r
y = angular_quadrature_points(k,2) * r
z = angular_quadrature_points(k,3) * r
accu += weights_angular_points(k) * (x**2 + y**2 + z**2)
enddo
print*, accu
end
subroutine routine
implicit none
integer :: i
double precision :: accu(2)
accu = 0.d0
do i = 1, N_det
call debug_det(psi_det(1,1,i),N_int)
enddo
do i = 1, nucl_num
accu(1) += integral_density_alpha_knowles_becke_per_atom(i)
accu(2) += integral_density_beta_knowles_becke_per_atom(i)
@ -19,6 +61,18 @@ subroutine routine
print*,'Nalpha = ',elec_alpha_num
print*,'accu(2) = ',accu(2)
print*,'Nalpha = ',elec_beta_num
accu = 0.d0
do i = 1, mo_tot_num
accu(1) += one_body_dm_mo_alpha_average(i,i)
accu(2) += one_body_dm_mo_beta_average(i,i)
enddo
print*,' '
print*,' '
print*,'accu(1) = ',accu(1)
print*,'accu(2) = ',accu(2)
end

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

15
plugins/MRPT/MRPT_Utils.main.irp.f
View File

@ -46,19 +46,6 @@ end
subroutine routine_2
implicit none
integer :: i
do i = 1, n_core_inact_orb
print*,fock_core_inactive_total(i,1,1),fock_core_inactive(i)
enddo
double precision :: accu
accu = 0.d0
do i = 1, n_act_orb
integer :: j_act_orb
j_act_orb = list_act(i)
accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb,1)
print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb,1),one_body_dm_mo_beta(j_act_orb,j_act_orb,1)
enddo
print*,'accu = ',accu
provide electronic_psi_ref_average_value
end

46
plugins/MRPT_Utils/MRMP2_density.irp.f Normal file
View File

@ -0,0 +1,46 @@
BEGIN_PROVIDER [double precision, MRMP2_density, (mo_tot_num_align, mo_tot_num)]
implicit none
integer :: i,j,k,l
double precision :: accu, mp2_dm(mo_tot_num)
MRMP2_density = one_body_dm_mo
call give_2h2p_density(mp2_dm)
accu = 0.d0
do i = 1, n_virt_orb
j = list_virt(i)
accu += mp2_dm(j)
MRMP2_density(j,j)+= mp2_dm(j)
enddo
END_PROVIDER
subroutine give_2h2p_density(mp2_density_diag_alpha_beta)
implicit none
double precision, intent(out) :: mp2_density_diag_alpha_beta(mo_tot_num)
integer :: i,j,k,l,m
integer :: iorb,jorb,korb,lorb
double precision :: get_mo_bielec_integral
double precision :: direct_int
double precision :: coef_double
mp2_density_diag_alpha_beta = 0.d0
do k = 1, n_virt_orb
korb = list_virt(k)
do i = 1, n_inact_orb
iorb = list_inact(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
do l = 1, n_virt_orb
lorb = list_virt(l)
direct_int = get_mo_bielec_integral(iorb,jorb,korb,lorb ,mo_integrals_map)
coef_double = direct_int/(fock_core_inactive_total_spin_trace(iorb,1) + fock_core_inactive_total_spin_trace(jorb,1) &
-fock_virt_total_spin_trace(korb,1) - fock_virt_total_spin_trace(lorb,1))
mp2_density_diag_alpha_beta(korb) += coef_double * coef_double
enddo
enddo
enddo
print*, mp2_density_diag_alpha_beta(korb)
enddo
end

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

@ -293,7 +293,8 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
END_PROVIDER
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
&BEGIN_PROVIDER [ double precision, two_anhil_one_creat_norm, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
implicit none
integer :: i,j
integer :: ispin,jspin,kspin
@ -344,6 +345,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) = norm_out(state_target)
enddo
enddo
enddo
@ -355,7 +357,54 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
BEGIN_PROVIDER [ double precision, two_anhil_one_creat_spin_average, (n_act_orb,n_act_orb,n_act_orb,N_states)]
implicit none
integer :: i,j
integer :: ispin,jspin,kspin
integer :: orb_i, hole_particle_i,spin_exc_i
integer :: orb_j, hole_particle_j,spin_exc_j
integer :: orb_k, hole_particle_k,spin_exc_k
double precision :: norm_out(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states)
double precision :: accu
do iorb = 1,n_act_orb
orb_i = list_act(iorb)
do jorb = 1, n_act_orb
orb_j = list_act(jorb)
do korb = 1, n_act_orb
orb_k = list_act(korb)
do state_target = 1, N_states
accu = 0.d0
do ispin = 1,2
do jspin = 1,2
do kspin = 1,2
two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) += two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target)* &
two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target)
accu += two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target)
enddo
enddo
enddo
two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) = two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) /accu
enddo
enddo
enddo
enddo
deallocate(psi_in_out,psi_in_out_coef)
END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
&BEGIN_PROVIDER [ double precision, two_creat_one_anhil_norm, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
implicit none
integer :: i,j
integer :: ispin,jspin,kspin
@ -406,6 +455,8 @@ implicit none
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) = norm_out(state_target)
! print*, norm_out(state_target)
enddo
enddo
enddo
@ -415,6 +466,51 @@ implicit none
enddo
deallocate(psi_in_out,psi_in_out_coef)
END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creat_one_anhil_spin_average, (n_act_orb,n_act_orb,n_act_orb,N_states)]
implicit none
integer :: i,j
integer :: ispin,jspin,kspin
integer :: orb_i, hole_particle_i,spin_exc_i
integer :: orb_j, hole_particle_j,spin_exc_j
integer :: orb_k, hole_particle_k,spin_exc_k
double precision :: norm_out(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states),accu
do iorb = 1,n_act_orb
orb_i = list_act(iorb)
do jorb = 1, n_act_orb
orb_j = list_act(jorb)
do korb = 1, n_act_orb
orb_k = list_act(korb)
do state_target = 1, N_states
accu = 0.d0
do ispin = 1,2
do jspin = 1,2
do kspin = 1,2
two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) += two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) * &
two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target)
accu += two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target)
print*, accu
enddo
enddo
enddo
two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) = two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) / accu
enddo
enddo
enddo
enddo
deallocate(psi_in_out,psi_in_out_coef)
END_PROVIDER
!BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,N_states)]
@ -1071,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

42
plugins/MRPT_Utils/ezfio_interface.irp.f
View File

@ -1,42 +0,0 @@
! DO NOT MODIFY BY HAND
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
! from file /home/giner/qp_bis/quantum_package/src/MRPT_Utils/EZFIO.cfg
BEGIN_PROVIDER [ logical, do_third_order_1h1p ]
implicit none
BEGIN_DOC
! If true, compute the third order contribution for the 1h1p
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_mrpt_utils_do_third_order_1h1p(has)
if (has) then
call ezfio_get_mrpt_utils_do_third_order_1h1p(do_third_order_1h1p)
else
print *, 'mrpt_utils/do_third_order_1h1p not found in EZFIO file'
stop 1
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, save_heff_eigenvectors ]
implicit none
BEGIN_DOC
! If true, save the eigenvectors of the dressed matrix at the end of the MRPT calculation
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_mrpt_utils_save_heff_eigenvectors(has)
if (has) then
call ezfio_get_mrpt_utils_save_heff_eigenvectors(save_heff_eigenvectors)
else
print *, 'mrpt_utils/save_heff_eigenvectors not found in EZFIO file'
stop 1
endif
END_PROVIDER

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

@ -121,14 +121,19 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
delta_e(i_state) = 1.d+20
enddo
else
call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e)
if(degree_scalar.eq.1)then
delta_e = 1.d+20
endif
! print*, 'delta_e',delta_e
!!!!!!!!!!!!! SHIFTED BK
! double precision :: hjj
! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
! delta_e(1) = CI_electronic_energy(1) - hjj
! delta_e(1) = electronic_psi_ref_average_value(1) - hjj
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
endif
hij_array(index_i) = hialpha
! print*, 'hialpha ',hialpha
do i_state = 1,N_states
delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
enddo

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

@ -34,43 +34,44 @@
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
enddo
second_order_pt_new_1h(i_state) = accu(i_state)
enddo
print*, '1h = ',accu
! 1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1p(delta_ij_tmp,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1p(i_state) = accu(i_state)
enddo
print*, '1p = ',accu
!! 1p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_1p(delta_ij_tmp,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_1p(i_state) = accu(i_state)
!enddo
!print*, '1p = ',accu
! 1h1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1h1p(delta_ij_tmp,N_det)
double precision :: e_corr_from_1h1p_singles(N_states)
!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles)
!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
enddo
print*, '1h1p = ',accu
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_1h1p(delta_ij_tmp,N_det)
!double precision :: e_corr_from_1h1p_singles(N_states)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_1h1p(i_state) = accu(i_state)
!enddo
!print*, '1h1p = ',accu
! 1h1p third order
if(do_third_order_1h1p)then
@ -83,75 +84,80 @@
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
enddo
print*, '1h1p(3)',accu
endif
! 2h
delta_ij_tmp = 0.d0
call H_apply_mrpt_2h(delta_ij_tmp,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2h(i_state) = accu(i_state)
enddo
print*, '2h = ',accu
!! 2h
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2h(delta_ij_tmp,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_2h(i_state) = accu(i_state)
!enddo
!print*, '2h = ',accu
! 2p
delta_ij_tmp = 0.d0
call H_apply_mrpt_2p(delta_ij_tmp,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2p(i_state) = accu(i_state)
enddo
print*, '2p = ',accu
!! 2p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2p(delta_ij_tmp,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_2p(i_state) = accu(i_state)
!enddo
!print*, '2p = ',accu
! 1h2p
delta_ij_tmp = 0.d0
!call give_1h2p_contrib(delta_ij_tmp)
call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1h2p(i_state) = accu(i_state)
enddo
print*, '1h2p = ',accu
!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_1h2p(i_state) = accu(i_state)
!enddo
!print*, '1h2p = ',accu
! 2h1p
delta_ij_tmp = 0.d0
!! 2h1p
!delta_ij_tmp = 0.d0
!call give_2h1p_contrib(delta_ij_tmp)
call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2h1p(i_state) = accu(i_state)
enddo
print*, '2h1p = ',accu
!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:)
!enddo
!second_order_pt_new_2h1p(i_state) = accu(i_state)
!enddo
!print*, '2h1p = ',accu
! 2h2p
!! 2h2p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det)
!accu = 0.d0
@ -178,10 +184,13 @@
! total
print*, ''
print*, 'total dressing'
print*, ''
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
do j = i_state, N_det
accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
enddo
@ -223,7 +232,7 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ]
BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ]
BEGIN_DOC
@ -245,7 +254,7 @@ END_PROVIDER
integer, allocatable :: iorder(:)
! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors
do j=1,min(N_states_diag,N_det)
do j=1,min(N_states,N_det)
do i=1,N_det
CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j)
enddo
@ -267,7 +276,7 @@ END_PROVIDER
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
Hmatrix_dressed_pt2_new_symmetrized,size(H_matrix_all_dets,1),N_det)
CI_electronic_energy(:) = 0.d0
if (s2_eig) then
i_state = 0
@ -276,8 +285,10 @@ END_PROVIDER
good_state_array = .False.
call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
N_det,size(eigenvectors,1))
print*,'N_det',N_det
do j=1,N_det
! Select at least n_states states with S^2 values closed to "expected_s2"
print*, s2_eigvalues(j),expected_s2
if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
i_state +=1
index_good_state_array(i_state) = j
@ -291,10 +302,10 @@ END_PROVIDER
! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
enddo
CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(index_good_state_array(j))
CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
enddo
i_other_state = 0
do j = 1, N_det
@ -304,10 +315,10 @@ END_PROVIDER
exit
endif
do i=1,N_det
CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
enddo
CI_electronic_energy(i_state+i_other_state) = eigenvalues(j)
CI_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
CI_dressed_pt2_new_electronic_energy(i_state+i_other_state) = eigenvalues(j)
CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
enddo
else
@ -322,10 +333,10 @@ END_PROVIDER
print*,''
do j=1,min(N_states_diag,N_det)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy(j) = eigenvalues(j)
CI_eigenvectors_s2(j) = s2_eigvalues(j)
CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j)
CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j)
enddo
endif
deallocate(index_good_state_array,good_state_array)
@ -336,9 +347,9 @@ END_PROVIDER
! Select the "N_states_diag" states of lowest energy
do j=1,min(N_det,N_states_diag)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy(j) = eigenvalues(j)
CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j)
enddo
endif
deallocate(eigenvectors,eigenvalues)
@ -358,7 +369,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
character*(8) :: st
call write_time(output_determinants)
do j=1,N_states_diag
CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
CI_dressed_pt2_new_energy(j) = CI_dressed_pt2_new_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))
call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))

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

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

@ -354,7 +354,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
kspin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
do i_state = 1, N_states
delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state)
! delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state)
delta_e_act(i_state) += two_anhil_one_creat_spin_average(i_particle_act,i_hole_act,j_hole_act,i_state)
enddo
else if (n_holes_act == 1 .and. n_particles_act == 2) then
@ -369,7 +370,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
j_particle_act = particle_list_practical(2,2)
do i_state = 1, N_states
delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state)
! delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state)
delta_e_act(i_state) += 0.5d0 * (two_creat_one_anhil_spin_average(i_particle_act,j_particle_act,i_hole_act,i_state) &
+two_creat_one_anhil_spin_average(j_particle_act,i_particle_act,i_hole_act,i_state))
enddo
else if (n_holes_act == 3 .and. n_particles_act == 0) then

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

14
plugins/Psiref_CAS/psi_ref.irp.f
View File

@ -72,10 +72,20 @@ END_PROVIDER
&BEGIN_PROVIDER [double precision, psi_ref_average_value, (N_states)]
implicit none
integer :: i,j
call u_0_H_u_0(electronic_psi_ref_average_value,psi_ref_coef,N_det_ref,psi_ref,N_int,N_states,psi_det_size)
electronic_psi_ref_average_value = psi_energy
do i = 1, N_states
psi_ref_average_value(i) = electronic_psi_ref_average_value(i) + nuclear_repulsion
psi_ref_average_value(i) = psi_energy(i) + nuclear_repulsion
enddo
double precision :: accu,hij
accu = 0.d0
do i = 1, N_det_ref
do j = 1, N_det_ref
call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
accu += psi_ref_coef(i,1) * psi_ref_coef(j,1) * hij
enddo
enddo
electronic_psi_ref_average_value(1) = accu
psi_ref_average_value(1) = electronic_psi_ref_average_value(1) + nuclear_repulsion
END_PROVIDER
BEGIN_PROVIDER [double precision, norm_psi_ref, (N_states)]

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

45
src/AO_Basis/ao_overlap.irp.f
View File

@ -129,3 +129,48 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_inv, (ao_num_align, ao_num) ]
implicit none
BEGIN_DOC
! Inverse of the overlap matrix
END_DOC
call invert_matrix(ao_overlap, size(ao_overlap,1), ao_num, ao_overlap_inv, size(ao_overlap_inv,1))
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_overlap_inv_1_2, (ao_num_align,ao_num)]
implicit none
integer :: i,j,k,l
double precision :: eigvalues(ao_num),eigvectors(ao_num_align, ao_num)
call lapack_diag(eigvalues,eigvectors,ao_overlap,ao_num_align,ao_num)
ao_overlap_inv_1_2 = 0.d0
double precision :: a_n
do i = 1, ao_num
a_n = 1.d0/dsqrt(eigvalues(i))
if(a_n.le.1.d-10)cycle
do j = 1, ao_num
do k = 1, ao_num
ao_overlap_inv_1_2(k,j) += eigvectors(k,i) * eigvectors(j,i) * a_n
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_overlap_1_2, (ao_num_align,ao_num)]
implicit none
integer :: i,j,k,l
double precision :: eigvalues(ao_num),eigvectors(ao_num_align, ao_num)
call lapack_diag(eigvalues,eigvectors,ao_overlap,ao_num_align,ao_num)
ao_overlap_1_2 = 0.d0
double precision :: a_n
do i = 1, ao_num
a_n = dsqrt(eigvalues(i))
do j = 1, ao_num
do k = 1, ao_num
ao_overlap_1_2(k,j) += eigvectors(k,i) * eigvectors(j,i) * a_n
enddo
enddo
enddo
END_PROVIDER

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

41
src/Determinants/print_wf.irp.f
View File

@ -32,29 +32,28 @@ subroutine routine
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'phase = ',phase
! if(degree == 1)then
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! if(s1 == 1)then
! norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
! else
! norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
! endif
if(degree == 1)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
if(s1 == 1)then
norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
! double precision :: hmono,hdouble
! call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
! print*,'hmono = ',hmono
! print*,'hdouble = ',hdouble
! print*,'hmono+hdouble = ',hmono+hdouble
! print*,'hij = ',hij
! else
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! print*,'s2',s2
! print*,'h2,p2 = ',h2,p2
double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
print*,'hmono = ',hmono
print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij
else if (degree == 2)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
! endif
endif
print*,'<Ref| H |D_I> = ',hij
endif
print*,'amplitude = ',psi_coef(i,1)/psi_coef(1,1)

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

92
src/MO_Basis/cholesky_mo.irp.f
View File

@ -50,12 +50,88 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
deallocate(W,work)
end
!subroutine svd_mo(n,m,P,LDP,C,LDC)
!implicit none
!BEGIN_DOC
! Singular value decomposition of the AO Density matrix
!
! n : Number of AOs
! m : Number of MOs
!
! P(LDP,n) : Density matrix in AO basis
!
! C(LDC,m) : MOs
!
! tol_in : tolerance
!
! rank : Nomber of local MOs (output)
!
!END_DOC
!integer, intent(in) :: n,m, LDC, LDP
!double precision, intent(in) :: P(LDP,n)
!double precision, intent(out) :: C(LDC,m)
!integer :: info
!integer :: i,k
!integer :: ipiv(n)
!double precision:: tol
!double precision, allocatable :: W(:,:), work(:)
!allocate(W(LDC,n),work(2*n))
!call svd(P,LDP,C,LDC,W,size(W,1),m,n)
!deallocate(W,work)
!end
subroutine svd_mo(n,m,P,LDP,C,LDC)
implicit none
BEGIN_DOC
! Singular value decomposition of the AO Density matrix
!
! n : Number of AOs
!
! m : Number of MOs
!
! P(LDP,n) : Density matrix in AO basis
!
! C(LDC,m) : MOs
!
END_DOC
integer, intent(in) :: n,m, LDC, LDP
double precision, intent(in) :: P(LDP,n)
double precision, intent(out) :: C(LDC,m)
integer :: info
integer :: i,k
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:), D(:)
allocate(W(LDC,n),work(2*n),D(n))
print*, ''
do i = 1, n
print*, P(i,i)
enddo
call svd(P,LDP,C,LDC,D,W,size(W,1),m,n)
double precision :: accu
accu = 0.d0
print*, 'm',m
do i = 1, m
print*, D(i)
accu += D(i)
enddo
print*,'Sum of D',accu
deallocate(W,work)
end
subroutine svd_mo_new(n,m,m_physical,P,LDP,C,LDC)
implicit none
BEGIN_DOC
! Singular value decomposition of the AO Density matrix
!
! n : Number of AOs
! m : Number of MOs
!
@ -68,7 +144,7 @@ subroutine svd_mo(n,m,P,LDP,C,LDC)
! rank : Nomber of local MOs (output)
!
END_DOC
integer, intent(in) :: n,m, LDC, LDP
integer, intent(in) :: n,m,m_physical, LDC, LDP
double precision, intent(in) :: P(LDP,n)
double precision, intent(out) :: C(LDC,m)
@ -76,10 +152,18 @@ subroutine svd_mo(n,m,P,LDP,C,LDC)
integer :: i,k
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:)
double precision, allocatable :: W(:,:), work(:), D(:)
allocate(W(LDC,n),work(2*n))
call svd(P,LDP,C,LDC,W,size(W,1),m,n)
allocate(W(LDC,n),work(2*n),D(n))
call svd(P,LDP,C,LDC,D,W,size(W,1),m_physical,n)
double precision :: accu
accu = 0.d0
print*, 'm',m_physical
do i = 1, m_physical
print*, D(i)
accu += D(i)
enddo
print*,'Sum of D',accu
deallocate(W,work)
end

122
src/MO_Basis/mos.irp.f
View File

@ -181,24 +181,146 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
! SC
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap,size(ao_overlap,1), &
mo_coef, size(mo_coef,1), &
0.d0, SC, ao_num_align)
! A.CS
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
SC, size(SC,1), &
0.d0, T, mo_tot_num_align)
! SC.A.CS
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC,size(SC,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
! C(S.A.S)C
! SC.A.CS
deallocate(T,SC)
end
subroutine mo_to_ao_s_inv_1_2(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the AO basis using the S^{-1} matrix
! S^{-1} C A C^{+} S^{-1}
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
double precision, allocatable :: T(:,:), SC_inv_1_2(:,:)
allocate ( SC_inv_1_2(ao_num_align,mo_tot_num) )
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
! SC_inv_1_2 = S^{-1}C
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap_inv_1_2,size(ao_overlap_inv_1_2,1), &
mo_coef, size(mo_coef,1), &
0.d0, SC_inv_1_2, ao_num_align)
! T = A.(SC_inv_1_2)^{+}
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
SC_inv_1_2, size(SC_inv_1_2,1), &
0.d0, T, mo_tot_num_align)
! SC_inv_1_2.A.CS
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC_inv_1_2,size(SC_inv_1_2,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
! C(S.A.S)C
! SC_inv_1_2.A.CS
deallocate(T,SC_inv_1_2)
end
subroutine mo_to_ao_s_1_2(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the AO basis using the S^{-1} matrix
! S^{-1} C A C^{+} S^{-1}
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
double precision, allocatable :: T(:,:), SC_1_2(:,:)
allocate ( SC_1_2(ao_num_align,mo_tot_num) )
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
! SC_1_2 = S^{-1}C
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap_1_2,size(ao_overlap_1_2,1), &
mo_coef, size(mo_coef,1), &
0.d0, SC_1_2, ao_num_align)
! T = A.(SC_1_2)^{+}
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
SC_1_2, size(SC_1_2,1), &
0.d0, T, mo_tot_num_align)
! SC_1_2.A.CS
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC_1_2,size(SC_1_2,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
! C(S.A.S)C
! SC_1_2.A.CS
deallocate(T,SC_1_2)
end
subroutine mo_to_ao_s_inv(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the AO basis using the S^{-1} matrix
! S^{-1} C A C^{+} S^{-1}
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
double precision, allocatable :: T(:,:), SC_inv(:,:)
allocate ( SC_inv(ao_num_align,mo_tot_num) )
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
! SC_inv = S^{-1}C
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap_inv,size(ao_overlap_inv,1), &
mo_coef, size(mo_coef,1), &
0.d0, SC_inv, ao_num_align)
! T = A.(SC_inv)^{+}
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
SC_inv, size(SC_inv,1), &
0.d0, T, mo_tot_num_align)
! SC_inv.A.CS
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC_inv,size(SC_inv,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
! C(S.A.S)C
! SC_inv.A.CS
deallocate(T,SC_inv)
end
subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC

16
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)') ''
@ -272,21 +272,13 @@ subroutine give_all_mos_at_r(r,mos_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
call give_specific_mos_at_r(r,mos_array, mo_coef)
end
subroutine give_specific_mos_at_r(r,mos_array, mo_coef_specific)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(in) :: mo_coef_specific(ao_num_align, mo_tot_num)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision :: aos_array(ao_num),accu
integer :: i,j
call give_all_aos_at_r(r,aos_array)
do i = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
accu += mo_coef_specific(j,i) * aos_array(j)
accu += mo_coef(j,i) * aos_array(j)
enddo
mos_array(i) = accu
enddo

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

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

@ -19,6 +19,10 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
double precision,allocatable :: A_tmp(:,:)
allocate (A_tmp(LDA,n))
print*, ''
do i = 1, n
print*, A(i,i)
enddo
A_tmp = A
! Find optimal size for temp arrays
@ -26,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))
@ -149,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 )
@ -293,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)

18
src/Utils/angular_integration.irp.f
View File

@ -4,7 +4,7 @@ BEGIN_PROVIDER [integer, degree_max_integration_lebedev]
! needed for the angular integration according to LEBEDEV formulae
END_DOC
implicit none
degree_max_integration_lebedev= 15
degree_max_integration_lebedev= 13
END_PROVIDER
@ -644,14 +644,14 @@ END_PROVIDER
weights_angular_integration_lebedev(16) = 0.016604069565742d0
weights_angular_integration_lebedev(17) = 0.016604069565742d0
weights_angular_integration_lebedev(18) = 0.016604069565742d0
weights_angular_integration_lebedev(19) = -0.029586038961039d0
weights_angular_integration_lebedev(20) = -0.029586038961039d0
weights_angular_integration_lebedev(21) = -0.029586038961039d0
weights_angular_integration_lebedev(22) = -0.029586038961039d0
weights_angular_integration_lebedev(23) = -0.029586038961039d0
weights_angular_integration_lebedev(24) = -0.029586038961039d0
weights_angular_integration_lebedev(25) = -0.029586038961039d0
weights_angular_integration_lebedev(26) = -0.029586038961039d0
weights_angular_integration_lebedev(19) = 0.029586038961039d0
weights_angular_integration_lebedev(20) = 0.029586038961039d0
weights_angular_integration_lebedev(21) = 0.029586038961039d0
weights_angular_integration_lebedev(22) = 0.029586038961039d0
weights_angular_integration_lebedev(23) = 0.029586038961039d0
weights_angular_integration_lebedev(24) = 0.029586038961039d0
weights_angular_integration_lebedev(25) = 0.029586038961039d0
weights_angular_integration_lebedev(26) = 0.029586038961039d0
weights_angular_integration_lebedev(27) = 0.026576207082159d0
weights_angular_integration_lebedev(28) = 0.026576207082159d0
weights_angular_integration_lebedev(29) = 0.026576207082159d0

5
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))
@ -9,3 +10,7 @@ double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
double precision, parameter :: thresh = 1.d-15
double precision, parameter :: cx_lda = -0.73855876638202234d0
double precision, parameter :: c_2_4_3 = 2.5198420997897464d0
double precision, parameter :: cst_lda = -0.93052573634909996d0
double precision, parameter :: c_4_3 = 1.3333333333333333d0

19
src/Utils/invert.irp.f Normal file
View File

@ -0,0 +1,19 @@
subroutine invert_matrix(A,LDA,na,A_inv,LDA_inv)
implicit none
double precision, intent(in) :: A (LDA,na)
integer, intent(in) :: LDA, LDA_inv
integer, intent(in) :: na
double precision, intent(out) :: A_inv (LDA_inv,na)
double precision :: work(LDA_inv*max(na,64))
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: work
integer :: inf
integer :: ipiv(LDA_inv)
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
integer :: lwork
A_inv(1:na,1:na) = A(1:na,1:na)
call dgetrf(na, na, A_inv, LDA_inv, ipiv, inf )
lwork = SIZE(work)
call dgetri(na, A_inv, LDA_inv, ipiv, work, lwork, inf )
end

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