mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Merge branch 'master' of https://github.com/scemama/quantum_package
Conflicts: config/ifort.cfg plugins/FOBOCI/H_apply_dressed_autonom.irp.f
This commit is contained in:
commit
67fd40626d
@ -31,13 +31,14 @@ OPENMP : 1 ; Append OpenMP flags
|
|||||||
# -ftz : Flushes denormal results to zero
|
# -ftz : Flushes denormal results to zero
|
||||||
#
|
#
|
||||||
[OPT]
|
[OPT]
|
||||||
FCFLAGS : -O2 -xHost -ip -ftz
|
FC : -traceback
|
||||||
|
FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g -traceback
|
||||||
|
|
||||||
# Profiling flags
|
# Profiling flags
|
||||||
#################
|
#################
|
||||||
#
|
#
|
||||||
[PROFILE]
|
[PROFILE]
|
||||||
FC : -p -g
|
FC : -p -g -traceback
|
||||||
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
|
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
|
||||||
|
|
||||||
# Debugging flags
|
# Debugging flags
|
||||||
@ -58,6 +59,6 @@ IRPF90_FLAGS : --openmp
|
|||||||
#################
|
#################
|
||||||
#
|
#
|
||||||
[OPENMP]
|
[OPENMP]
|
||||||
FC : -openmp
|
FC : -qopenmp
|
||||||
IRPF90_FLAGS : --openmp
|
IRPF90_FLAGS : --openmp
|
||||||
|
|
||||||
|
7
configure
vendored
7
configure
vendored
@ -46,7 +46,12 @@ if len(sys.argv) != 3:
|
|||||||
# \_| |_ (_) |_) (_| | | | | | (_)
|
# \_| |_ (_) |_) (_| | | | | | (_)
|
||||||
#
|
#
|
||||||
|
|
||||||
QP_ROOT = os.getcwd()
|
try:
|
||||||
|
QP_ROOT = os.environ["QP_ROOT"]
|
||||||
|
except KeyError:
|
||||||
|
QP_ROOT = os.getcwd()
|
||||||
|
os.environ["QP_ROOT"] = QP_ROOT
|
||||||
|
|
||||||
QP_ROOT_BIN = join(QP_ROOT, "bin")
|
QP_ROOT_BIN = join(QP_ROOT, "bin")
|
||||||
QP_ROOT_LIB = join(QP_ROOT, "lib")
|
QP_ROOT_LIB = join(QP_ROOT, "lib")
|
||||||
QP_ROOT_INSTALL = join(QP_ROOT, "install")
|
QP_ROOT_INSTALL = join(QP_ROOT, "install")
|
||||||
|
@ -49,20 +49,20 @@ let t_to_string = function
|
|||||||
| None -> assert false
|
| None -> assert false
|
||||||
;;
|
;;
|
||||||
|
|
||||||
let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename =
|
let set ~core ~inact ~act ~virt ~del =
|
||||||
|
|
||||||
Ezfio.set_file ezfio_filename ;
|
let mo_tot_num =
|
||||||
if not (Ezfio.has_mo_basis_mo_tot_num ()) then
|
Ezfio.get_mo_basis_mo_tot_num ()
|
||||||
failwith "mo_basis/mo_tot_num not found" ;
|
in
|
||||||
|
|
||||||
let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in
|
|
||||||
let n_int =
|
let n_int =
|
||||||
try N_int_number.of_int (Ezfio.get_determinants_n_int ())
|
try N_int_number.of_int (Ezfio.get_determinants_n_int ())
|
||||||
with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num
|
with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num
|
||||||
in
|
in
|
||||||
|
|
||||||
|
|
||||||
let mo_class = Array.init mo_tot_num ~f:(fun i -> None) in
|
let mo_class =
|
||||||
|
Array.init mo_tot_num ~f:(fun i -> None)
|
||||||
|
in
|
||||||
|
|
||||||
(* Check input data *)
|
(* Check input data *)
|
||||||
let apply_class l =
|
let apply_class l =
|
||||||
@ -196,6 +196,49 @@ let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_fi
|
|||||||
|> Ezfio.set_bitmasks_cas;
|
|> Ezfio.set_bitmasks_cas;
|
||||||
;;
|
;;
|
||||||
|
|
||||||
|
let get () =
|
||||||
|
|
||||||
|
let mo_tot_num =
|
||||||
|
Ezfio.get_mo_basis_mo_tot_num ()
|
||||||
|
in
|
||||||
|
let n_int =
|
||||||
|
try N_int_number.of_int (Ezfio.get_determinants_n_int ())
|
||||||
|
with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num
|
||||||
|
in
|
||||||
|
|
||||||
|
let bitmasks =
|
||||||
|
match Input.Bitmasks.read () with
|
||||||
|
| Some x -> x
|
||||||
|
| None -> failwith "No data to print"
|
||||||
|
in
|
||||||
|
assert (bitmasks.Input.Bitmasks.n_mask_gen |> Bitmask_number.to_int = 1);
|
||||||
|
assert (bitmasks.Input.Bitmasks.n_mask_cas |> Bitmask_number.to_int = 1);
|
||||||
|
|
||||||
|
let (generators,cas) =
|
||||||
|
Bitlist.of_int64_array bitmasks.Input.Bitmasks.generators,
|
||||||
|
Bitlist.of_int64_array bitmasks.Input.Bitmasks.cas
|
||||||
|
in
|
||||||
|
|
||||||
|
Printf.printf "MO : %d\n" mo_tot_num;
|
||||||
|
Printf.printf "n_int: %d\n" (N_int_number.to_int n_int);
|
||||||
|
Printf.printf "Gen : %s\nCAS : %s\n"
|
||||||
|
(Bitlist.to_string generators)
|
||||||
|
(Bitlist.to_string cas)
|
||||||
|
|
||||||
|
;;
|
||||||
|
|
||||||
|
let run ~print ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename =
|
||||||
|
|
||||||
|
Ezfio.set_file ezfio_filename ;
|
||||||
|
if not (Ezfio.has_mo_basis_mo_tot_num ()) then
|
||||||
|
failwith "mo_basis/mo_tot_num not found" ;
|
||||||
|
|
||||||
|
if print then
|
||||||
|
get ()
|
||||||
|
else
|
||||||
|
set ~core ~inact ~act ~virt ~del
|
||||||
|
;;
|
||||||
|
|
||||||
let ezfio_file =
|
let ezfio_file =
|
||||||
let failure filename =
|
let failure filename =
|
||||||
eprintf "'%s' is not an EZFIO file.\n%!" filename;
|
eprintf "'%s' is not an EZFIO file.\n%!" filename;
|
||||||
@ -240,6 +283,7 @@ let spec =
|
|||||||
+> flag "act" (optional string) ~doc:"range Range of active orbitals"
|
+> flag "act" (optional string) ~doc:"range Range of active orbitals"
|
||||||
+> flag "virt" (optional string) ~doc:"range Range of virtual orbitals"
|
+> flag "virt" (optional string) ~doc:"range Range of virtual orbitals"
|
||||||
+> flag "del" (optional string) ~doc:"range Range of deleted orbitals"
|
+> flag "del" (optional string) ~doc:"range Range of deleted orbitals"
|
||||||
|
+> flag "print" no_arg ~doc:" Print the current masks"
|
||||||
+> anon ("ezfio_filename" %: ezfio_file)
|
+> anon ("ezfio_filename" %: ezfio_file)
|
||||||
;;
|
;;
|
||||||
|
|
||||||
@ -251,7 +295,7 @@ let command =
|
|||||||
The range of MOs has the form : \"[36-53,72-107,126-131]\"
|
The range of MOs has the form : \"[36-53,72-107,126-131]\"
|
||||||
")
|
")
|
||||||
spec
|
spec
|
||||||
(fun core inact act virt del ezfio_filename () -> run ?core ?inact ?act ?virt ?del ezfio_filename )
|
(fun core inact act virt del print ezfio_filename () -> run ~print ?core ?inact ?act ?virt ?del ezfio_filename )
|
||||||
;;
|
;;
|
||||||
|
|
||||||
let () =
|
let () =
|
||||||
|
@ -273,7 +273,8 @@ subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,ipr
|
|||||||
integer,parameter :: size_max = 3072
|
integer,parameter :: size_max = 3072
|
||||||
|
|
||||||
integer, intent(in) :: Ndet_generators
|
integer, intent(in) :: Ndet_generators
|
||||||
double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
|
double precision, intent(in) :: E_ref
|
||||||
|
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
|
||||||
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
||||||
|
|
||||||
integer ,intent(in) :: i_generator
|
integer ,intent(in) :: i_generator
|
||||||
@ -437,8 +438,9 @@ subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_g
|
|||||||
|
|
||||||
|
|
||||||
integer, intent(in) :: Ndet_generators
|
integer, intent(in) :: Ndet_generators
|
||||||
|
double precision, intent(in) :: E_ref
|
||||||
|
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
|
||||||
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
||||||
double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
|
|
||||||
|
|
||||||
|
|
||||||
integer :: i_generator, nmax
|
integer :: i_generator, nmax
|
||||||
|
@ -2,20 +2,25 @@ use bitmasks
|
|||||||
BEGIN_SHELL [ /usr/bin/env python ]
|
BEGIN_SHELL [ /usr/bin/env python ]
|
||||||
from generate_h_apply import *
|
from generate_h_apply import *
|
||||||
|
|
||||||
s = H_apply_zmq("FCI")
|
s = H_apply("FCI")
|
||||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||||
s.unset_openmp()
|
#s.unset_openmp()
|
||||||
print s
|
print s
|
||||||
|
|
||||||
|
#s = H_apply("FCI_PT2")
|
||||||
|
#s.set_perturbation("epstein_nesbet_2x2")
|
||||||
|
#s.unset_openmp()
|
||||||
|
#print s
|
||||||
|
|
||||||
s = H_apply_zmq("FCI_PT2")
|
s = H_apply_zmq("FCI_PT2")
|
||||||
s.set_perturbation("epstein_nesbet_2x2")
|
s.set_perturbation("epstein_nesbet_2x2")
|
||||||
s.unset_openmp()
|
s.unset_openmp()
|
||||||
print s
|
print s
|
||||||
|
|
||||||
s = H_apply_zmq("FCI_no_skip")
|
s = H_apply("FCI_no_skip")
|
||||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||||
s.unset_skip()
|
s.unset_skip()
|
||||||
s.unset_openmp()
|
#s.unset_openmp()
|
||||||
print s
|
print s
|
||||||
|
|
||||||
s = H_apply("FCI_mono")
|
s = H_apply("FCI_mono")
|
||||||
|
@ -2,3 +2,16 @@
|
|||||||
type: double precision
|
type: double precision
|
||||||
doc: Calculated energy
|
doc: Calculated energy
|
||||||
interface: ezfio
|
interface: ezfio
|
||||||
|
|
||||||
|
[thresh_mrcc]
|
||||||
|
type: Threshold
|
||||||
|
doc: Threshold on the convergence of the MRCC energy
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 1.e-5
|
||||||
|
|
||||||
|
[n_it_mrcc_max]
|
||||||
|
type: Strictly_positive_int
|
||||||
|
doc: Maximum number of MRCC iterations
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 10
|
||||||
|
|
||||||
|
@ -1,13 +1,109 @@
|
|||||||
program mrcc
|
program mrcc
|
||||||
implicit none
|
implicit none
|
||||||
if (.not.read_wf) then
|
double precision, allocatable :: energy(:)
|
||||||
print *, 'read_wf has to be true.'
|
allocate (energy(N_states))
|
||||||
stop 1
|
|
||||||
endif
|
read_wf = .True.
|
||||||
|
SOFT_TOUCH read_wf
|
||||||
call print_cas_coefs
|
call print_cas_coefs
|
||||||
call run_mrcc
|
call set_generators_bitmasks_as_holes_and_particles
|
||||||
|
call run(N_states,energy)
|
||||||
|
if(do_pt2_end)then
|
||||||
|
call run_pt2(N_states,energy)
|
||||||
|
endif
|
||||||
|
deallocate(energy)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
subroutine run(N_st,energy)
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: N_st
|
||||||
|
double precision, intent(out) :: energy(N_st)
|
||||||
|
|
||||||
|
integer :: i
|
||||||
|
|
||||||
|
double precision :: E_new, E_old, delta_e
|
||||||
|
integer :: iteration
|
||||||
|
double precision :: E_past(4), lambda
|
||||||
|
E_new = 0.d0
|
||||||
|
delta_E = 1.d0
|
||||||
|
iteration = 0
|
||||||
|
lambda = 1.d0
|
||||||
|
do while (delta_E > thresh_mrcc)
|
||||||
|
iteration += 1
|
||||||
|
print *, '==========================='
|
||||||
|
print *, 'MRCC Iteration', iteration
|
||||||
|
print *, '==========================='
|
||||||
|
print *, ''
|
||||||
|
E_old = sum(ci_energy_dressed)
|
||||||
|
call write_double(6,ci_energy_dressed(1),"MRCC energy")
|
||||||
|
call diagonalize_ci_dressed(lambda)
|
||||||
|
E_new = sum(ci_energy_dressed)
|
||||||
|
delta_E = dabs(E_new - E_old)
|
||||||
|
call save_wavefunction
|
||||||
|
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
|
||||||
|
if (iteration > n_it_mrcc_max) then
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
|
||||||
|
energy(:) = ci_energy_dressed(:)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine run_pt2(N_st,energy)
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,k
|
||||||
|
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
|
||||||
|
integer, intent(in) :: N_st
|
||||||
|
double precision, intent(in) :: energy(N_st)
|
||||||
|
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
|
||||||
|
pt2 = 0.d0
|
||||||
|
|
||||||
|
print*,'Last iteration only to compute the PT2'
|
||||||
|
threshold_selectors = 1.d0
|
||||||
|
threshold_generators = 0.999d0
|
||||||
|
|
||||||
|
N_det_generators = lambda_mrcc_pt2(0) + N_det_cas
|
||||||
|
do i=1,N_det_cas
|
||||||
|
do k=1,N_int
|
||||||
|
psi_det_generators(k,1,i) = psi_ref(k,1,i)
|
||||||
|
psi_det_generators(k,2,i) = psi_ref(k,2,i)
|
||||||
|
enddo
|
||||||
|
do k=1,N_st
|
||||||
|
psi_coef_generators(i,k) = psi_ref_coef(i,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do i=N_det_cas+1,N_det_generators
|
||||||
|
j = lambda_mrcc_pt2(i)
|
||||||
|
do k=1,N_int
|
||||||
|
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
|
||||||
|
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
|
||||||
|
enddo
|
||||||
|
do k=1,N_st
|
||||||
|
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
|
||||||
|
|
||||||
|
|
||||||
|
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
|
||||||
|
print *, 'Final step'
|
||||||
|
print *, 'N_det = ', N_det
|
||||||
|
print *, 'N_states = ', N_states
|
||||||
|
print *, 'PT2 = ', pt2
|
||||||
|
print *, 'E = ', energy
|
||||||
|
print *, 'E+PT2 = ', energy+pt2
|
||||||
|
print *, '-----'
|
||||||
|
|
||||||
|
|
||||||
|
call ezfio_set_full_ci_energy_pt2(energy+pt2)
|
||||||
|
deallocate(pt2,norm_pert)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
subroutine print_cas_coefs
|
subroutine print_cas_coefs
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
@ -18,7 +114,7 @@ subroutine print_cas_coefs
|
|||||||
print *, psi_cas_coef(i,:)
|
print *, psi_cas_coef(i,:)
|
||||||
call debug_det(psi_cas(1,1,i),N_int)
|
call debug_det(psi_cas(1,1,i),N_int)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
call write_double(6,ci_energy(1),"Initial CI energy")
|
call write_double(6,ci_energy(1),"Initial CI energy")
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
91
plugins/MRCC_CASSD/mrcc_noiter.irp.f
Normal file
91
plugins/MRCC_CASSD/mrcc_noiter.irp.f
Normal file
@ -0,0 +1,91 @@
|
|||||||
|
program mrcc_noiter
|
||||||
|
implicit none
|
||||||
|
double precision, allocatable :: energy(:)
|
||||||
|
allocate (energy(N_states))
|
||||||
|
read_wf = .True.
|
||||||
|
threshold_generators = .9999d0
|
||||||
|
SOFT_TOUCH read_wf threshold_generators
|
||||||
|
call print_cas_coefs
|
||||||
|
call set_generators_bitmasks_as_holes_and_particles
|
||||||
|
call run(N_states,energy)
|
||||||
|
if(do_pt2_end)then
|
||||||
|
call run_pt2(N_states,energy)
|
||||||
|
endif
|
||||||
|
deallocate(energy)
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine run(N_st,energy)
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: N_st
|
||||||
|
double precision, intent(out) :: energy(N_st)
|
||||||
|
integer :: i,j
|
||||||
|
do j=1,N_states_diag
|
||||||
|
do i=1,N_det
|
||||||
|
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
SOFT_TOUCH psi_coef ci_energy_dressed
|
||||||
|
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
|
||||||
|
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
|
||||||
|
call save_wavefunction
|
||||||
|
energy(:) = ci_energy_dressed(:)
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine run_pt2(N_st,energy)
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,k
|
||||||
|
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
|
||||||
|
integer, intent(in) :: N_st
|
||||||
|
double precision, intent(in) :: energy(N_st)
|
||||||
|
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
|
||||||
|
pt2 = 0.d0
|
||||||
|
|
||||||
|
print*,'Last iteration only to compute the PT2'
|
||||||
|
threshold_selectors = 1.d0
|
||||||
|
threshold_generators = 0.999d0
|
||||||
|
|
||||||
|
N_det_generators = lambda_mrcc_pt2(0)
|
||||||
|
do i=1,N_det_generators
|
||||||
|
j = lambda_mrcc_pt2(i)
|
||||||
|
do k=1,N_int
|
||||||
|
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
|
||||||
|
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
|
||||||
|
enddo
|
||||||
|
do k=1,N_st
|
||||||
|
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
|
||||||
|
|
||||||
|
|
||||||
|
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
|
||||||
|
print *, 'Final step'
|
||||||
|
print *, 'N_det = ', N_det
|
||||||
|
print *, 'N_states = ', N_states
|
||||||
|
print *, 'PT2 = ', pt2
|
||||||
|
print *, 'E = ', energy
|
||||||
|
print *, 'E+PT2 = ', energy+pt2
|
||||||
|
print *, '-----'
|
||||||
|
|
||||||
|
call ezfio_set_full_ci_energy_pt2(energy+pt2)
|
||||||
|
deallocate(pt2,norm_pert)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine print_cas_coefs
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer :: i,j
|
||||||
|
print *, 'CAS'
|
||||||
|
print *, '==='
|
||||||
|
do i=1,N_det_cas
|
||||||
|
print *, psi_cas_coef(i,:)
|
||||||
|
call debug_det(psi_cas(1,1,i),N_int)
|
||||||
|
enddo
|
||||||
|
call write_double(6,ci_energy(1),"Initial CI energy")
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -3,19 +3,19 @@ BEGIN_SHELL [ /usr/bin/env python ]
|
|||||||
from generate_h_apply import *
|
from generate_h_apply import *
|
||||||
|
|
||||||
s = H_apply("mrcc")
|
s = H_apply("mrcc")
|
||||||
s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
|
s.data["parameters"] = ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
|
||||||
s.data["declarations"] += """
|
s.data["declarations"] += """
|
||||||
integer, intent(in) :: Ndet_ref,Ndet_non_ref
|
integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
|
||||||
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, Ndet_ref)
|
||||||
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
|
double precision, intent(in) :: delta_ii_(Nstates, Ndet_ref)
|
||||||
"""
|
"""
|
||||||
s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
|
s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
|
||||||
s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
|
s.data["params_post"] += ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
|
||||||
s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
|
s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
|
||||||
s.data["decls_main"] += """
|
s.data["decls_main"] += """
|
||||||
integer, intent(in) :: Ndet_ref,Ndet_non_ref
|
integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates
|
||||||
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
|
||||||
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
|
double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref)
|
||||||
"""
|
"""
|
||||||
s.data["finalization"] = ""
|
s.data["finalization"] = ""
|
||||||
s.data["copy_buffer"] = ""
|
s.data["copy_buffer"] = ""
|
||||||
@ -24,27 +24,18 @@ s.data["size_max"] = "3072"
|
|||||||
print s
|
print s
|
||||||
|
|
||||||
|
|
||||||
s = H_apply("mrcepa")
|
|
||||||
s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
|
|
||||||
s.data["declarations"] += """
|
|
||||||
integer, intent(in) :: Ndet_ref,Ndet_non_ref
|
|
||||||
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
|
||||||
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
|
|
||||||
"""
|
|
||||||
s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
|
|
||||||
s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
|
|
||||||
s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
|
|
||||||
s.data["decls_main"] += """
|
|
||||||
integer, intent(in) :: Ndet_ref,Ndet_non_ref
|
|
||||||
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
|
||||||
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
|
|
||||||
"""
|
|
||||||
s.data["finalization"] = ""
|
|
||||||
s.data["copy_buffer"] = ""
|
|
||||||
s.data["generate_psi_guess"] = ""
|
|
||||||
s.data["size_max"] = "3072"
|
|
||||||
# print s
|
|
||||||
|
|
||||||
|
s = H_apply("mrcc_PT2")
|
||||||
|
s.energy = "ci_electronic_energy_dressed"
|
||||||
|
s.set_perturbation("epstein_nesbet_2x2")
|
||||||
|
s.unset_openmp()
|
||||||
|
print s
|
||||||
|
|
||||||
|
#s = H_apply_zmq("mrcc_PT2")
|
||||||
|
#s.energy = "ci_electronic_energy_dressed"
|
||||||
|
#s.set_perturbation("epstein_nesbet_2x2")
|
||||||
|
#s.unset_openmp()
|
||||||
|
#print s
|
||||||
|
|
||||||
END_SHELL
|
END_SHELL
|
||||||
|
|
||||||
|
@ -189,10 +189,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
|
|||||||
! Davidson iterations
|
! Davidson iterations
|
||||||
! ===================
|
! ===================
|
||||||
|
|
||||||
|
|
||||||
|
integer :: iteration
|
||||||
converged = .False.
|
converged = .False.
|
||||||
|
|
||||||
do while (.not.converged)
|
do while (.not.converged)
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
!$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st)
|
!$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st)
|
||||||
do k=1,N_st
|
do k=1,N_st
|
||||||
@ -206,6 +206,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
|
|||||||
|
|
||||||
do iter=1,davidson_sze_max-1
|
do iter=1,davidson_sze_max-1
|
||||||
|
|
||||||
|
|
||||||
! Compute W_k = H |u_k>
|
! Compute W_k = H |u_k>
|
||||||
! ----------------------
|
! ----------------------
|
||||||
|
|
||||||
|
@ -14,17 +14,17 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
|
|||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
|
subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
integer, intent(in) :: i_generator,n_selected, Nint, iproc
|
integer, intent(in) :: i_generator,n_selected, Nint, iproc
|
||||||
integer, intent(in) :: Ndet_ref, Ndet_non_ref
|
integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
|
||||||
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
|
||||||
double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
|
double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
|
||||||
|
|
||||||
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
||||||
integer :: i,j,k,l
|
integer :: i,j,k,l,m
|
||||||
integer :: degree_alpha(psi_det_size)
|
integer :: degree_alpha(psi_det_size)
|
||||||
integer :: idx_alpha(0:psi_det_size)
|
integer :: idx_alpha(0:psi_det_size)
|
||||||
logical :: good, fullMatch
|
logical :: good, fullMatch
|
||||||
@ -32,10 +32,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
|
|||||||
integer(bit_kind) :: tq(Nint,2,n_selected)
|
integer(bit_kind) :: tq(Nint,2,n_selected)
|
||||||
integer :: N_tq, c_ref ,degree
|
integer :: N_tq, c_ref ,degree
|
||||||
|
|
||||||
double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
|
double precision :: hIk, hla, hIl, dIk(Nstates), dka(Nstates), dIa(Nstates)
|
||||||
double precision, allocatable :: dIa_hla(:,:)
|
double precision, allocatable :: dIa_hla(:,:)
|
||||||
double precision :: haj, phase, phase2
|
double precision :: haj, phase, phase2
|
||||||
double precision :: f(N_states), ci_inv(N_states)
|
double precision :: f(Nstates), ci_inv(Nstates)
|
||||||
integer :: exc(0:2,2,2)
|
integer :: exc(0:2,2,2)
|
||||||
integer :: h1,h2,p1,p2,s1,s2
|
integer :: h1,h2,p1,p2,s1,s2
|
||||||
integer(bit_kind) :: tmp_det(Nint,2)
|
integer(bit_kind) :: tmp_det(Nint,2)
|
||||||
@ -46,10 +46,16 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
|
|||||||
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
|
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
|
||||||
integer,allocatable :: idx_miniList(:)
|
integer,allocatable :: idx_miniList(:)
|
||||||
integer :: N_miniList, ni, leng
|
integer :: N_miniList, ni, leng
|
||||||
|
double precision, allocatable :: hij_cache(:)
|
||||||
|
|
||||||
|
integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
|
||||||
|
integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
|
||||||
|
integer :: mobiles(2), smallerlist
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
leng = max(N_det_generators, N_det_non_ref)
|
leng = max(N_det_generators, N_det_non_ref)
|
||||||
allocate(miniList(Nint, 2, leng), idx_miniList(leng))
|
allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref))
|
||||||
|
|
||||||
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
|
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
|
||||||
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
|
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
|
||||||
@ -58,27 +64,129 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
|
|||||||
return
|
return
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
allocate(ptr_microlist(0:mo_tot_num*2+1), &
|
||||||
|
N_microlist(0:mo_tot_num*2) )
|
||||||
|
allocate( microlist(Nint,2,N_minilist*4), &
|
||||||
|
idx_microlist(N_minilist*4))
|
||||||
|
|
||||||
|
if(key_mask(1,1) /= 0) then
|
||||||
|
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
|
||||||
|
call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
|
||||||
|
else
|
||||||
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
|
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
|
||||||
|
end if
|
||||||
|
|
||||||
allocate (dIa_hla(N_states,Ndet_non_ref))
|
|
||||||
|
|
||||||
|
deallocate(microlist, idx_microlist)
|
||||||
|
|
||||||
|
allocate (dIa_hla(Nstates,Ndet_non_ref))
|
||||||
|
|
||||||
! |I>
|
! |I>
|
||||||
|
|
||||||
! |alpha>
|
! |alpha>
|
||||||
|
|
||||||
if(N_tq > 0) then
|
if(N_tq > 0) then
|
||||||
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
|
call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint)
|
||||||
|
if(N_minilist == 0) return
|
||||||
|
|
||||||
|
|
||||||
|
if(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!!
|
||||||
|
allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist))
|
||||||
|
|
||||||
|
allocate( microlist(Nint,2,N_minilist*4), &
|
||||||
|
idx_microlist(N_minilist*4))
|
||||||
|
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
|
||||||
|
|
||||||
|
|
||||||
|
do i=0,mo_tot_num*2
|
||||||
|
do k=ptr_microlist(i),ptr_microlist(i+1)-1
|
||||||
|
idx_microlist(k) = idx_minilist(idx_microlist(k))
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do l=1,N_microlist(0)
|
||||||
|
do k=1,Nint
|
||||||
|
microlist_zero(k,1,l) = microlist(k,1,l)
|
||||||
|
microlist_zero(k,2,l) = microlist(k,2,l)
|
||||||
|
enddo
|
||||||
|
idx_microlist_zero(l) = idx_microlist(l)
|
||||||
|
enddo
|
||||||
|
end if
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
|
||||||
do i_alpha=1,N_tq
|
|
||||||
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
|
|
||||||
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
|
|
||||||
|
|
||||||
|
|
||||||
|
do i_alpha=1,N_tq
|
||||||
|
if(key_mask(1,1) /= 0) then
|
||||||
|
call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint)
|
||||||
|
|
||||||
|
if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
|
||||||
|
smallerlist = mobiles(1)
|
||||||
|
else
|
||||||
|
smallerlist = mobiles(2)
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
do l=0,N_microlist(smallerlist)-1
|
||||||
|
microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l)
|
||||||
|
idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l)
|
||||||
|
end do
|
||||||
|
|
||||||
|
call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha)
|
||||||
|
do j=1,idx_alpha(0)
|
||||||
|
idx_alpha(j) = idx_microlist_zero(idx_alpha(j))
|
||||||
|
end do
|
||||||
|
|
||||||
|
|
||||||
|
! i = 1
|
||||||
|
! j = 2
|
||||||
|
! do j = 2, idx_alpha_tmp(0)
|
||||||
|
! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
! m = j
|
||||||
|
!
|
||||||
|
! idx_alpha(0) = idx_alpha_tmp(0)
|
||||||
|
!
|
||||||
|
! do l = 1, idx_alpha(0)
|
||||||
|
! if(j > idx_alpha_tmp(0)) then
|
||||||
|
! k = i
|
||||||
|
! i += 1
|
||||||
|
! else if(i >= m) then
|
||||||
|
! k = j
|
||||||
|
! j += 1
|
||||||
|
! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then
|
||||||
|
! k = i
|
||||||
|
! i += 1
|
||||||
|
! else
|
||||||
|
! k = j
|
||||||
|
! j += 1
|
||||||
|
! end if
|
||||||
|
! ! k=l
|
||||||
|
! idx_alpha(l) = idx_alpha_tmp(k)
|
||||||
|
! degree_alpha(l) = degree_alpha_tmp(k)
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
else
|
||||||
|
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
|
||||||
do j=1,idx_alpha(0)
|
do j=1,idx_alpha(0)
|
||||||
idx_alpha(j) = idx_miniList(idx_alpha(j))
|
idx_alpha(j) = idx_miniList(idx_alpha(j))
|
||||||
end do
|
end do
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
|
||||||
|
! do j=1,idx_alpha(0)
|
||||||
|
! idx_alpha(j) = idx_miniList(idx_alpha(j))
|
||||||
|
! end do
|
||||||
|
!print *, idx_alpha(:idx_alpha(0))
|
||||||
|
|
||||||
|
do l_sd=1,idx_alpha(0)
|
||||||
|
k_sd = idx_alpha(l_sd)
|
||||||
|
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
|
||||||
|
enddo
|
||||||
|
|
||||||
! |I>
|
! |I>
|
||||||
do i_I=1,N_det_ref
|
do i_I=1,N_det_ref
|
||||||
@ -88,20 +196,36 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
|
|||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
|
||||||
do i_state=1,N_states
|
do i_state=1,Nstates
|
||||||
dIa(i_state) = 0.d0
|
dIa(i_state) = 0.d0
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! <I| <> |alpha>
|
! <I| <> |alpha>
|
||||||
do k_sd=1,idx_alpha(0)
|
do k_sd=1,idx_alpha(0)
|
||||||
|
|
||||||
|
! Loop if lambda == 0
|
||||||
|
logical :: loop
|
||||||
|
loop = .True.
|
||||||
|
do i_state=1,Nstates
|
||||||
|
if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then
|
||||||
|
loop = .False.
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
if (loop) then
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
|
||||||
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
|
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
|
||||||
if (degree > 2) then
|
if (degree > 2) then
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
|
||||||
! <I| /k\ |alpha>
|
! <I| /k\ |alpha>
|
||||||
! <I|H|k>
|
! <I|H|k>
|
||||||
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
|
hIk = hij_mrcc(idx_alpha(k_sd),i_I)
|
||||||
do i_state=1,N_states
|
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
|
||||||
|
do i_state=1,Nstates
|
||||||
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
|
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
|
||||||
enddo
|
enddo
|
||||||
! |l> = Exc(k -> alpha) |I>
|
! |l> = Exc(k -> alpha) |I>
|
||||||
@ -111,87 +235,81 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
|
|||||||
tmp_det(k,1) = psi_ref(k,1,i_I)
|
tmp_det(k,1) = psi_ref(k,1,i_I)
|
||||||
tmp_det(k,2) = psi_ref(k,2,i_I)
|
tmp_det(k,2) = psi_ref(k,2,i_I)
|
||||||
enddo
|
enddo
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
|
|
||||||
|
|
||||||
! Particle
|
logical :: ok
|
||||||
iint = ishft(p1-1,-bit_kind_shift) + 1
|
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
|
||||||
ipos = p1-ishft((iint-1),bit_kind_shift)-1
|
if(.not. ok) cycle
|
||||||
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
|
|
||||||
if (degree_alpha(k_sd) == 2) then
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
|
|
||||||
endif
|
|
||||||
|
|
||||||
! <I| \l/ |alpha>
|
! <I| \l/ |alpha>
|
||||||
do i_state=1,N_states
|
do i_state=1,Nstates
|
||||||
dka(i_state) = 0.d0
|
dka(i_state) = 0.d0
|
||||||
enddo
|
enddo
|
||||||
do l_sd=k_sd+1,idx_alpha(0)
|
do l_sd=k_sd+1,idx_alpha(0)
|
||||||
|
|
||||||
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
|
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
|
||||||
if (degree == 0) then
|
if (degree == 0) then
|
||||||
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
|
|
||||||
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
|
loop = .True.
|
||||||
do i_state=1,N_states
|
do i_state=1,Nstates
|
||||||
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
|
if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then
|
||||||
enddo
|
loop = .False.
|
||||||
exit
|
exit
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do i_state=1,N_states
|
if (.not.loop) then
|
||||||
|
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
|
||||||
|
hIl = hij_mrcc(idx_alpha(l_sd),i_I)
|
||||||
|
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
|
||||||
|
do i_state=1,Nstates
|
||||||
|
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do i_state=1,Nstates
|
||||||
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
|
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do i_state=1,N_states
|
do i_state=1,Nstates
|
||||||
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
|
ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state)
|
||||||
enddo
|
enddo
|
||||||
do l_sd=1,idx_alpha(0)
|
do l_sd=1,idx_alpha(0)
|
||||||
k_sd = idx_alpha(l_sd)
|
k_sd = idx_alpha(l_sd)
|
||||||
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
|
hla = hij_cache(k_sd)
|
||||||
do i_state=1,N_states
|
! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
|
||||||
|
do i_state=1,Nstates
|
||||||
dIa_hla(i_state,k_sd) = dIa(i_state) * hla
|
dIa_hla(i_state,k_sd) = dIa(i_state) * hla
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
call omp_set_lock( psi_ref_lock(i_I) )
|
call omp_set_lock( psi_ref_lock(i_I) )
|
||||||
|
do i_state=1,Nstates
|
||||||
|
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
|
||||||
do l_sd=1,idx_alpha(0)
|
do l_sd=1,idx_alpha(0)
|
||||||
k_sd = idx_alpha(l_sd)
|
k_sd = idx_alpha(l_sd)
|
||||||
do i_state=1,N_states
|
delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
|
||||||
delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd)
|
delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
|
||||||
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
|
|
||||||
delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state)
|
|
||||||
else
|
|
||||||
delta_ii_(i_I,i_state) = 0.d0
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
|
else
|
||||||
|
delta_ii_(i_state,i_I) = 0.d0
|
||||||
|
do l_sd=1,idx_alpha(0)
|
||||||
|
k_sd = idx_alpha(l_sd)
|
||||||
|
delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
enddo
|
enddo
|
||||||
call omp_unset_lock( psi_ref_lock(i_I) )
|
call omp_unset_lock( psi_ref_lock(i_I) )
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
deallocate (dIa_hla)
|
!deallocate (dIa_hla,hij_cache)
|
||||||
deallocate(miniList, idx_miniList)
|
!deallocate(miniList, idx_miniList)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ]
|
|
||||||
&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
|
|
||||||
&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
|
|
||||||
&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
|
|
||||||
gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
|
|
||||||
gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
|
|
||||||
call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
|
|
||||||
call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
|
subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
|
||||||
@ -224,6 +342,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
|
|||||||
N_tq = 0
|
N_tq = 0
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
i_loop : do i=1,N_selected
|
i_loop : do i=1,N_selected
|
||||||
if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then
|
if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then
|
||||||
cycle
|
cycle
|
||||||
@ -253,8 +372,84 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
|
|||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
|
||||||
|
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: i_generator,n_selected, Nint
|
||||||
|
|
||||||
|
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
||||||
|
integer :: i,j,k,m
|
||||||
|
logical :: is_in_wavefunction
|
||||||
|
integer :: degree(psi_det_size)
|
||||||
|
integer :: idx(0:psi_det_size)
|
||||||
|
logical :: good
|
||||||
|
|
||||||
|
integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
|
||||||
|
integer, intent(out) :: N_tq
|
||||||
|
|
||||||
|
|
||||||
|
integer :: nt,ni
|
||||||
|
logical, external :: is_connected_to
|
||||||
|
|
||||||
|
|
||||||
|
integer(bit_kind),intent(in) :: microlist(Nint,2,*)
|
||||||
|
integer,intent(in) :: ptr_microlist(0:*)
|
||||||
|
integer,intent(in) :: N_microlist(0:*)
|
||||||
|
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
|
||||||
|
|
||||||
|
integer :: mobiles(2), smallerlist
|
||||||
|
|
||||||
|
N_tq = 0
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
i_loop : do i=1,N_selected
|
||||||
|
call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint)
|
||||||
|
if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
|
||||||
|
smallerlist = mobiles(1)
|
||||||
|
else
|
||||||
|
smallerlist = mobiles(2)
|
||||||
|
end if
|
||||||
|
|
||||||
|
if(N_microlist(smallerlist) > 0) then
|
||||||
|
if(is_connected_to(det_buffer(1,1,i), microlist(1,1,ptr_microlist(smallerlist)), Nint, N_microlist(smallerlist))) then
|
||||||
|
cycle
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
if(N_microlist(0) > 0) then
|
||||||
|
if(is_connected_to(det_buffer(1,1,i), microlist, Nint, N_microlist(0))) then
|
||||||
|
cycle
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Select determinants that are triple or quadruple excitations
|
||||||
|
! from the ref
|
||||||
|
good = .True.
|
||||||
|
call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
|
||||||
|
!good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
|
||||||
|
do k=1,idx(0)
|
||||||
|
if (degree(k) < 3) then
|
||||||
|
good = .False.
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
if (good) then
|
||||||
|
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
|
||||||
|
N_tq += 1
|
||||||
|
do k=1,N_int
|
||||||
|
tq(k,1,N_tq) = det_buffer(k,1,i)
|
||||||
|
tq(k,2,N_tq) = det_buffer(k,2,i)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
enddo i_loop
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
@ -11,12 +11,13 @@ subroutine mrcc_iterations
|
|||||||
|
|
||||||
double precision :: E_new, E_old, delta_e
|
double precision :: E_new, E_old, delta_e
|
||||||
integer :: iteration,i_oscillations
|
integer :: iteration,i_oscillations
|
||||||
double precision :: E_past(4)
|
double precision :: E_past(4), lambda
|
||||||
E_new = 0.d0
|
E_new = 0.d0
|
||||||
delta_E = 1.d0
|
delta_E = 1.d0
|
||||||
iteration = 0
|
iteration = 0
|
||||||
j = 1
|
j = 1
|
||||||
i_oscillations = 0
|
i_oscillations = 0
|
||||||
|
lambda = 1.d0
|
||||||
do while (delta_E > 1.d-7)
|
do while (delta_E > 1.d-7)
|
||||||
iteration += 1
|
iteration += 1
|
||||||
print *, '==========================='
|
print *, '==========================='
|
||||||
@ -25,29 +26,18 @@ subroutine mrcc_iterations
|
|||||||
print *, ''
|
print *, ''
|
||||||
E_old = sum(ci_energy_dressed)
|
E_old = sum(ci_energy_dressed)
|
||||||
call write_double(6,ci_energy_dressed(1),"MRCC energy")
|
call write_double(6,ci_energy_dressed(1),"MRCC energy")
|
||||||
call diagonalize_ci_dressed
|
call diagonalize_ci_dressed(lambda)
|
||||||
E_new = sum(ci_energy_dressed)
|
E_new = sum(ci_energy_dressed)
|
||||||
delta_E = dabs(E_new - E_old)
|
delta_E = dabs(E_new - E_old)
|
||||||
|
! if (E_new > E_old) then
|
||||||
|
! lambda = lambda * 0.7d0
|
||||||
|
! else
|
||||||
|
! lambda = min(1.d0, lambda * 1.1d0)
|
||||||
|
! endif
|
||||||
|
! print *, 'energy lambda ', lambda
|
||||||
E_past(j) = E_new
|
E_past(j) = E_new
|
||||||
j +=1
|
j +=1
|
||||||
if(j>4)then
|
|
||||||
j=1
|
|
||||||
endif
|
|
||||||
if(iteration > 4) then
|
|
||||||
if(delta_E > 1.d-10)then
|
|
||||||
if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then
|
|
||||||
print*,'OSCILLATIONS !!!'
|
|
||||||
oscillations = .True.
|
|
||||||
i_oscillations +=1
|
|
||||||
lambda_mrcc_tmp = lambda_mrcc
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
call save_wavefunction
|
call save_wavefunction
|
||||||
! if (i_oscillations > 5) then
|
|
||||||
! exit
|
|
||||||
! endif
|
|
||||||
if (iteration > 200) then
|
if (iteration > 200) then
|
||||||
exit
|
exit
|
||||||
endif
|
endif
|
||||||
|
@ -1,115 +1,71 @@
|
|||||||
BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ]
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
|
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
|
||||||
&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ]
|
&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
|
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
|
||||||
END_DOC
|
END_DOC
|
||||||
integer :: i,k,j
|
integer :: i,k
|
||||||
double precision :: ihpsi(N_states), hii,delta_e_eff,ihpsi_current(N_states),hij
|
double precision :: ihpsi_current(N_states)
|
||||||
integer :: i_ok,i_pert,i_pert_count
|
integer :: i_pert_count
|
||||||
i_ok = 0
|
double precision :: hii, lambda_pert
|
||||||
|
integer :: N_lambda_mrcc_pt2
|
||||||
|
|
||||||
double precision :: phase_restart(N_states),tmp
|
|
||||||
do k = 1, N_states
|
|
||||||
phase_restart(k) = dsign(1.d0,psi_ref_coef_restart(1,k)/psi_ref_coef(1,k))
|
|
||||||
enddo
|
|
||||||
i_pert_count = 0
|
i_pert_count = 0
|
||||||
|
lambda_mrcc = 0.d0
|
||||||
|
N_lambda_mrcc_pt2 = 0
|
||||||
|
lambda_mrcc_pt2(0) = 0
|
||||||
|
|
||||||
do i=1,N_det_non_ref
|
do i=1,N_det_non_ref
|
||||||
call i_h_psi(psi_non_ref(1,1,i), psi_ref_restart, psi_ref_coef_restart, N_int, N_det_ref,&
|
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
|
||||||
size(psi_ref_coef_restart,1), n_states, ihpsi)
|
size(psi_ref_coef,1), N_states,ihpsi_current)
|
||||||
call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
|
call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
|
||||||
! TODO --- Test perturbatif ------
|
|
||||||
do k=1,N_states
|
do k=1,N_states
|
||||||
lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
|
|
||||||
! TODO : i_h_psi peut sortir de la boucle?
|
|
||||||
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
|
|
||||||
if (ihpsi_current(k) == 0.d0) then
|
if (ihpsi_current(k) == 0.d0) then
|
||||||
ihpsi_current(k) = 1.d-32
|
ihpsi_current(k) = 1.d-32
|
||||||
endif
|
endif
|
||||||
tmp = psi_non_ref_coef(i,k)/ihpsi_current(k)
|
lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) )
|
||||||
i_pert = 0
|
lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
|
||||||
! Perturbation only if 1st order < 0.5 x second order
|
if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then
|
||||||
if((ihpsi(k) * lambda_pert(k,i)) < 0.5d0 * psi_non_ref_coef_restart(i,k) )then
|
i_pert_count += 1
|
||||||
i_pert = 1
|
lambda_mrcc(k,i) = 0.d0
|
||||||
else
|
if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then
|
||||||
do j = 1, N_det_ref
|
N_lambda_mrcc_pt2 += 1
|
||||||
call i_H_j(psi_non_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
|
lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i
|
||||||
! Perturbation diverges when hij*tmp > 0.5
|
endif
|
||||||
if(dabs(hij * tmp).ge.0.5d0)then
|
|
||||||
i_pert_count +=1
|
|
||||||
i_pert = 1
|
|
||||||
exit
|
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
endif
|
|
||||||
if( i_pert == 1)then
|
|
||||||
pert_determinants(k,i) = i_pert
|
|
||||||
endif
|
|
||||||
if(pert_determinants(k,i) == 1)then
|
|
||||||
i_ok +=1
|
|
||||||
lambda_mrcc(k,i) = lambda_pert(k,i)
|
|
||||||
else
|
|
||||||
lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
! TODO --- Fin test perturbatif ------
|
lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2
|
||||||
enddo
|
|
||||||
!if(oscillations)then
|
|
||||||
! print*,'AVERAGING the lambda_mrcc with those of the previous iterations'
|
|
||||||
! do i = 1, N_det_non_ref
|
|
||||||
! do k = 1, N_states
|
|
||||||
|
|
||||||
! double precision :: tmp
|
|
||||||
! tmp = lambda_mrcc(k,i)
|
|
||||||
! lambda_mrcc(k,i) += lambda_mrcc_tmp(k,i)
|
|
||||||
! lambda_mrcc(k,i) = lambda_mrcc(k,i) * 0.5d0
|
|
||||||
! if(dabs(tmp - lambda_mrcc(k,i)).ge.1.d-9)then
|
|
||||||
! print*,''
|
|
||||||
! print*,'i = ',i
|
|
||||||
! print*,'psi_non_ref_coef(i,k) = ',psi_non_ref_coef(i,k)
|
|
||||||
! print*,'lambda_mrcc(k,i) = ',lambda_mrcc(k,i)
|
|
||||||
! print*,' tmp = ',tmp
|
|
||||||
! endif
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!endif
|
|
||||||
print*,'N_det_non_ref = ',N_det_non_ref
|
print*,'N_det_non_ref = ',N_det_non_ref
|
||||||
print*,'Number of Perturbatively treated determinants = ',i_ok
|
print*,'Number of ignored determinants = ',i_pert_count
|
||||||
print*,'i_pert_count = ',i_pert_count
|
|
||||||
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
|
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
|
||||||
|
print*,'lambda max = ',maxval(dabs(lambda_mrcc))
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, lambda_mrcc_tmp, (N_states,psi_det_size) ]
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
|
||||||
implicit none
|
implicit none
|
||||||
lambda_mrcc_tmp = 0.d0
|
BEGIN_DOC
|
||||||
|
! < ref | H | Non-ref > matrix
|
||||||
|
END_DOC
|
||||||
|
integer :: i_I, k_sd
|
||||||
|
do i_I=1,N_det_ref
|
||||||
|
do k_sd=1,N_det_non_ref
|
||||||
|
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,k_sd),N_int,hij_mrcc(k_sd,i_I))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ logical, oscillations ]
|
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
|
||||||
implicit none
|
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states,N_det_ref) ]
|
||||||
oscillations = .False.
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
|
|
||||||
!implicit none
|
|
||||||
!BEGIN_DOC
|
|
||||||
!! Dressing matrix in SD basis
|
|
||||||
!END_DOC
|
|
||||||
!delta_ij_non_ref = 0.d0
|
|
||||||
!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
|
|
||||||
!END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
|
||||||
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
|
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Dressing matrix in N_det basis
|
! Dressing matrix in N_det basis
|
||||||
@ -117,32 +73,7 @@ END_PROVIDER
|
|||||||
integer :: i,j,m
|
integer :: i,j,m
|
||||||
delta_ij = 0.d0
|
delta_ij = 0.d0
|
||||||
delta_ii = 0.d0
|
delta_ii = 0.d0
|
||||||
call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref)
|
call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref)
|
||||||
double precision :: max_delta
|
|
||||||
double precision :: accu
|
|
||||||
integer :: imax,jmax
|
|
||||||
max_delta = 0.d0
|
|
||||||
accu = 0.d0
|
|
||||||
do i = 1, N_det_ref
|
|
||||||
do j = 1, N_det_non_ref
|
|
||||||
accu += psi_non_ref_coef(j,1) * psi_ref_coef(i,1) * delta_ij(i,j,1)
|
|
||||||
if(dabs(delta_ij(i,j,1)).gt.max_delta)then
|
|
||||||
max_delta = dabs(delta_ij(i,j,1))
|
|
||||||
imax = i
|
|
||||||
jmax = j
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
print*,''
|
|
||||||
print*,''
|
|
||||||
print*,'<psi| Delta H |psi> = ',accu
|
|
||||||
print*,'MAX VAL OF DRESING = ',delta_ij(imax,jmax,1)
|
|
||||||
print*,'imax,jmax = ',imax,jmax
|
|
||||||
print*,'psi_ref_coef(imax,1) = ',psi_ref_coef(imax,1)
|
|
||||||
print*,'psi_non_ref_coef(jmax,1) = ',psi_non_ref_coef(jmax,1)
|
|
||||||
do i = 1, N_det_ref
|
|
||||||
print*,'delta_ii(i,1) = ',delta_ii(i,1)
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
|
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
|
||||||
@ -159,11 +90,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
|
|||||||
enddo
|
enddo
|
||||||
do ii = 1, N_det_ref
|
do ii = 1, N_det_ref
|
||||||
i =idx_ref(ii)
|
i =idx_ref(ii)
|
||||||
h_matrix_dressed(i,i,istate) += delta_ii(ii,istate)
|
h_matrix_dressed(i,i,istate) += delta_ii(istate,ii)
|
||||||
do jj = 1, N_det_non_ref
|
do jj = 1, N_det_non_ref
|
||||||
j =idx_non_ref(jj)
|
j =idx_non_ref(jj)
|
||||||
h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate)
|
h_matrix_dressed(i,j,istate) += delta_ij(istate,jj,ii)
|
||||||
h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate)
|
h_matrix_dressed(j,i,istate) += delta_ij(istate,jj,ii)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -252,18 +183,21 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
|
|||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
subroutine diagonalize_CI_dressed
|
subroutine diagonalize_CI_dressed(lambda)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Replace the coefficients of the CI states by the coefficients of the
|
! Replace the coefficients of the CI states by the coefficients of the
|
||||||
! eigenstates of the CI matrix
|
! eigenstates of the CI matrix
|
||||||
END_DOC
|
END_DOC
|
||||||
|
double precision, intent(in) :: lambda
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
do j=1,N_states_diag
|
do j=1,N_states_diag
|
||||||
do i=1,N_det
|
do i=1,N_det
|
||||||
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
|
psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j)
|
||||||
enddo
|
enddo
|
||||||
|
call normalize(psi_coef(1,j), N_det)
|
||||||
enddo
|
enddo
|
||||||
SOFT_TOUCH psi_coef
|
SOFT_TOUCH psi_coef
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -1,260 +0,0 @@
|
|||||||
use omp_lib
|
|
||||||
use bitmasks
|
|
||||||
|
|
||||||
subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
|
|
||||||
use bitmasks
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer, intent(in) :: i_generator,n_selected, Nint, iproc
|
|
||||||
integer, intent(in) :: Ndet_ref, Ndet_non_ref
|
|
||||||
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
|
|
||||||
double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
|
|
||||||
|
|
||||||
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
|
||||||
integer :: i,j,k,l
|
|
||||||
integer :: degree_alpha(psi_det_size)
|
|
||||||
integer :: idx_alpha(0:psi_det_size)
|
|
||||||
logical :: good, fullMatch
|
|
||||||
|
|
||||||
integer(bit_kind) :: tq(Nint,2,n_selected)
|
|
||||||
integer :: N_tq, c_ref ,degree
|
|
||||||
|
|
||||||
double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
|
|
||||||
double precision, allocatable :: dIa_hia(:,:)
|
|
||||||
double precision :: haj, phase, phase2
|
|
||||||
double precision :: f(N_states), ci_inv(N_states)
|
|
||||||
integer :: exc(0:2,2,2)
|
|
||||||
integer :: h1,h2,p1,p2,s1,s2
|
|
||||||
integer(bit_kind) :: tmp_det(Nint,2)
|
|
||||||
integer(bit_kind) :: tmp_det_0(Nint,2)
|
|
||||||
integer :: iint, ipos
|
|
||||||
integer :: i_state, i_sd, k_sd, l_sd, i_I, i_alpha
|
|
||||||
|
|
||||||
integer(bit_kind),allocatable :: miniList(:,:,:)
|
|
||||||
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
|
|
||||||
integer,allocatable :: idx_miniList(:)
|
|
||||||
integer :: N_miniList, ni, leng
|
|
||||||
integer(bit_kind) :: isum
|
|
||||||
|
|
||||||
double precision :: hia
|
|
||||||
integer, allocatable :: index_sorted(:)
|
|
||||||
|
|
||||||
|
|
||||||
leng = max(N_det_generators, N_det_non_ref)
|
|
||||||
allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det))
|
|
||||||
|
|
||||||
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
|
|
||||||
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
|
|
||||||
|
|
||||||
if(fullMatch) then
|
|
||||||
return
|
|
||||||
end if
|
|
||||||
|
|
||||||
|
|
||||||
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
|
|
||||||
|
|
||||||
allocate (dIa_hia(N_states,Ndet_non_ref))
|
|
||||||
|
|
||||||
! |I>
|
|
||||||
|
|
||||||
! |alpha>
|
|
||||||
|
|
||||||
if(N_tq > 0) then
|
|
||||||
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
|
|
||||||
end if
|
|
||||||
|
|
||||||
|
|
||||||
do i_alpha=1,N_tq
|
|
||||||
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
|
|
||||||
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
|
|
||||||
|
|
||||||
integer, external :: get_index_in_psi_det_sorted_bit
|
|
||||||
index_sorted = huge(-1)
|
|
||||||
do j=1,idx_alpha(0)
|
|
||||||
idx_alpha(j) = idx_miniList(idx_alpha(j))
|
|
||||||
index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j)
|
|
||||||
end do
|
|
||||||
|
|
||||||
! |I>
|
|
||||||
do i_I=1,N_det_ref
|
|
||||||
! Find triples and quadruple grand parents
|
|
||||||
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
|
|
||||||
if (degree > 4) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
|
|
||||||
do i_state=1,N_states
|
|
||||||
dIa(i_state) = 0.d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!TODO: MR
|
|
||||||
do i_sd=1,idx_alpha(0)
|
|
||||||
call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint)
|
|
||||||
if (degree > 2) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
|
|
||||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
||||||
tmp_det_0 = 0_bit_kind
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
|
|
||||||
if (degree == 2) then
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
|
|
||||||
endif
|
|
||||||
|
|
||||||
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia)
|
|
||||||
|
|
||||||
! <I| <> |alpha>
|
|
||||||
do k_sd=1,idx_alpha(0)
|
|
||||||
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
|
|
||||||
if (degree > 2) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
|
|
||||||
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint)
|
|
||||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
||||||
tmp_det = 0_bit_kind
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
|
|
||||||
if (degree == 2) then
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
|
|
||||||
endif
|
|
||||||
|
|
||||||
isum = 0_bit_kind
|
|
||||||
do iint = 1,N_int
|
|
||||||
isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) &
|
|
||||||
+ iand(tmp_det(iint,2), tmp_det_0(iint,2))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
if (isum /= 0_bit_kind) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
|
|
||||||
! <I| /k\ |alpha>
|
|
||||||
! <I|H|k>
|
|
||||||
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
|
|
||||||
do i_state=1,N_states
|
|
||||||
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
|
|
||||||
enddo
|
|
||||||
! |l> = Exc(k -> alpha) |I>
|
|
||||||
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
|
|
||||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
||||||
do k=1,N_int
|
|
||||||
tmp_det(k,1) = psi_ref(k,1,i_I)
|
|
||||||
tmp_det(k,2) = psi_ref(k,2,i_I)
|
|
||||||
enddo
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p1-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p1-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
|
|
||||||
if (degree == 2) then
|
|
||||||
! Hole (see list_to_bitstring)
|
|
||||||
iint = ishft(h2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = h2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
|
|
||||||
|
|
||||||
! Particle
|
|
||||||
iint = ishft(p2-1,-bit_kind_shift) + 1
|
|
||||||
ipos = p2-ishft((iint-1),bit_kind_shift)-1
|
|
||||||
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
|
|
||||||
endif
|
|
||||||
|
|
||||||
! <I| \l/ |alpha>
|
|
||||||
do i_state=1,N_states
|
|
||||||
dka(i_state) = 0.d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) )
|
|
||||||
! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint)
|
|
||||||
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl)
|
|
||||||
! do i_state=1,N_states
|
|
||||||
! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2
|
|
||||||
! enddo
|
|
||||||
|
|
||||||
do l_sd=1,idx_alpha(0)
|
|
||||||
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
|
|
||||||
if (degree == 0) then
|
|
||||||
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
|
|
||||||
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
|
|
||||||
do i_state=1,N_states
|
|
||||||
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
|
|
||||||
enddo
|
|
||||||
exit
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
do i_state=1,N_states
|
|
||||||
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do i_state=1,N_states
|
|
||||||
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
k_sd = idx_alpha(i_sd)
|
|
||||||
do i_state=1,N_states
|
|
||||||
dIa_hia(i_state,k_sd) = dIa(i_state) * hia
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call omp_set_lock( psi_ref_lock(i_I) )
|
|
||||||
do i_state=1,N_states
|
|
||||||
delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd)
|
|
||||||
|
|
||||||
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
|
|
||||||
delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state)
|
|
||||||
else
|
|
||||||
delta_ii_(i_I,i_state) = 0.d0
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
call omp_unset_lock( psi_ref_lock(i_I) )
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
enddo
|
|
||||||
deallocate (dIa_hia,index_sorted)
|
|
||||||
deallocate(miniList, idx_miniList)
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -1,97 +0,0 @@
|
|||||||
subroutine run_mrcepa
|
|
||||||
implicit none
|
|
||||||
call set_generators_bitmasks_as_holes_and_particles
|
|
||||||
call mrcepa_iterations
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine mrcepa_iterations
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i,j
|
|
||||||
|
|
||||||
double precision :: E_new, E_old, delta_e
|
|
||||||
integer :: iteration,i_oscillations
|
|
||||||
double precision :: E_past(4)
|
|
||||||
E_new = 0.d0
|
|
||||||
delta_E = 1.d0
|
|
||||||
iteration = 0
|
|
||||||
j = 1
|
|
||||||
i_oscillations = 0
|
|
||||||
do while (delta_E > 1.d-7)
|
|
||||||
iteration += 1
|
|
||||||
print *, '==========================='
|
|
||||||
print *, 'MRCEPA Iteration', iteration
|
|
||||||
print *, '==========================='
|
|
||||||
print *, ''
|
|
||||||
E_old = sum(ci_energy_dressed)
|
|
||||||
call write_double(6,ci_energy_dressed(1),"MRCEPA energy")
|
|
||||||
call diagonalize_ci_dressed
|
|
||||||
E_new = sum(ci_energy_dressed)
|
|
||||||
delta_E = dabs(E_new - E_old)
|
|
||||||
|
|
||||||
E_past(j) = E_new
|
|
||||||
j +=1
|
|
||||||
if(j>4)then
|
|
||||||
j=1
|
|
||||||
endif
|
|
||||||
if(iteration > 4) then
|
|
||||||
if(delta_E > 1.d-10)then
|
|
||||||
if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then
|
|
||||||
print*,'OSCILLATIONS !!!'
|
|
||||||
oscillations = .True.
|
|
||||||
i_oscillations +=1
|
|
||||||
lambda_mrcc_tmp = lambda_mrcc
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
call save_wavefunction
|
|
||||||
! if (i_oscillations > 5) then
|
|
||||||
! exit
|
|
||||||
! endif
|
|
||||||
if (iteration > 200) then
|
|
||||||
exit
|
|
||||||
endif
|
|
||||||
print*,'------------'
|
|
||||||
print*,'VECTOR'
|
|
||||||
do i = 1, N_det_ref
|
|
||||||
print*,''
|
|
||||||
print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
|
|
||||||
print*,'delta_ii(i,1) = ',delta_ii(i,1)
|
|
||||||
enddo
|
|
||||||
print*,'------------'
|
|
||||||
enddo
|
|
||||||
call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy")
|
|
||||||
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
|
|
||||||
call save_wavefunction
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine set_generators_bitmasks_as_holes_and_particles
|
|
||||||
implicit none
|
|
||||||
integer :: i,k
|
|
||||||
do k = 1, N_generators_bitmask
|
|
||||||
do i = 1, N_int
|
|
||||||
! Pure single part
|
|
||||||
generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha
|
|
||||||
generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha
|
|
||||||
generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta
|
|
||||||
generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta
|
|
||||||
|
|
||||||
! Double excitation
|
|
||||||
generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha
|
|
||||||
generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha
|
|
||||||
generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta
|
|
||||||
generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta
|
|
||||||
|
|
||||||
generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha
|
|
||||||
generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha
|
|
||||||
generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta
|
|
||||||
generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
touch generators_bitmask
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
end
|
|
@ -3,7 +3,7 @@ import perturbation
|
|||||||
END_SHELL
|
END_SHELL
|
||||||
|
|
||||||
|
|
||||||
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
|
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
|
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
|
||||||
@ -14,6 +14,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
|
|||||||
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
|
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
|
||||||
integer(bit_kind),intent(in) :: key_mask(Nint,2)
|
integer(bit_kind),intent(in) :: key_mask(Nint,2)
|
||||||
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
|
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
|
||||||
|
double precision, intent(in) :: electronic_energy(N_st)
|
||||||
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
|
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
|
||||||
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
|
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
|
||||||
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
|
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
|
||||||
@ -151,7 +152,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
|
|||||||
idx_microlist_zero(ptr_microlist(1)+l) = idx_microlist(ptr_microlist(smallerlist)+l)
|
idx_microlist_zero(ptr_microlist(1)+l) = idx_microlist(ptr_microlist(smallerlist)+l)
|
||||||
enddo
|
enddo
|
||||||
end if
|
end if
|
||||||
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
||||||
c_pert,e_2_pert,H_pert_diag,Nint,N_microlist(smallerlist)+N_microlist(0), &
|
c_pert,e_2_pert,H_pert_diag,Nint,N_microlist(smallerlist)+N_microlist(0), &
|
||||||
n_st,microlist_zero,idx_microlist_zero,N_microlist(smallerlist)+N_microlist(0))
|
n_st,microlist_zero,idx_microlist_zero,N_microlist(smallerlist)+N_microlist(0))
|
||||||
else
|
else
|
||||||
@ -160,11 +161,11 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
|
|||||||
cycle
|
cycle
|
||||||
end if
|
end if
|
||||||
|
|
||||||
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
||||||
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
! call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
||||||
! c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
! c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
||||||
|
|
||||||
do k = 1,N_st
|
do k = 1,N_st
|
||||||
@ -182,7 +183,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
|
|||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
|
subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
|
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
|
||||||
@ -193,6 +194,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
|
|||||||
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
|
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
|
||||||
integer(bit_kind),intent(in) :: key_mask(Nint,2)
|
integer(bit_kind),intent(in) :: key_mask(Nint,2)
|
||||||
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
|
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
|
||||||
|
double precision, intent(in) :: electronic_energy(N_st)
|
||||||
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
|
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
|
||||||
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
|
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
|
||||||
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
|
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
|
||||||
@ -241,7 +243,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
|
|||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
|
||||||
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
|
||||||
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
|
||||||
|
|
||||||
do k = 1,N_st
|
do k = 1,N_st
|
||||||
|
@ -29,11 +29,11 @@ subroutine pt2_epstein_nesbet ($arguments)
|
|||||||
|
|
||||||
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
||||||
do i =1,N_st
|
do i =1,N_st
|
||||||
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
|
if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then
|
||||||
c_pert(i) = -1.d0
|
c_pert(i) = -1.d0
|
||||||
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
|
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
|
||||||
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
|
else if (dabs(electronic_energy(i) - h) > 1.d-6) then
|
||||||
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
|
c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h)
|
||||||
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
|
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
|
||||||
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
||||||
else
|
else
|
||||||
@ -66,7 +66,6 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
|
|||||||
double precision :: i_H_psi_array(N_st)
|
double precision :: i_H_psi_array(N_st)
|
||||||
ASSERT (Nint == N_int)
|
ASSERT (Nint == N_int)
|
||||||
ASSERT (Nint > 0)
|
ASSERT (Nint > 0)
|
||||||
PROVIDE CI_electronic_energy
|
|
||||||
|
|
||||||
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
|
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
|
||||||
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
|
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
|
||||||
@ -74,7 +73,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
|
|||||||
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
||||||
do i =1,N_st
|
do i =1,N_st
|
||||||
if (i_H_psi_array(i) /= 0.d0) then
|
if (i_H_psi_array(i) /= 0.d0) then
|
||||||
delta_e = h - CI_electronic_energy(i)
|
delta_e = h - electronic_energy(i)
|
||||||
if (delta_e > 0.d0) then
|
if (delta_e > 0.d0) then
|
||||||
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
|
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
|
||||||
else
|
else
|
||||||
@ -167,7 +166,7 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
|
|||||||
!
|
!
|
||||||
! that can be repeated by repeating all the double excitations
|
! that can be repeated by repeating all the double excitations
|
||||||
!
|
!
|
||||||
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
|
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
|
||||||
!
|
!
|
||||||
! that could be repeated to this determinant.
|
! that could be repeated to this determinant.
|
||||||
!
|
!
|
||||||
@ -197,16 +196,16 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
|
|||||||
enddo
|
enddo
|
||||||
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
||||||
h = h + accu_e_corr
|
h = h + accu_e_corr
|
||||||
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
|
delta_E = 1.d0/(electronic_energy(1) - h)
|
||||||
|
|
||||||
|
|
||||||
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
|
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
|
||||||
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
|
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
|
||||||
|
|
||||||
do i =2,N_st
|
do i =2,N_st
|
||||||
H_pert_diag(i) = h
|
H_pert_diag(i) = h
|
||||||
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
|
if (dabs(electronic_energy(i) - h) > 1.d-6) then
|
||||||
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
|
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
|
||||||
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
|
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
|
||||||
else
|
else
|
||||||
c_pert(i) = i_H_psi_array(i)
|
c_pert(i) = i_H_psi_array(i)
|
||||||
@ -250,7 +249,7 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
|
|||||||
!
|
!
|
||||||
! that can be repeated by repeating all the double excitations
|
! that can be repeated by repeating all the double excitations
|
||||||
!
|
!
|
||||||
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
|
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
|
||||||
!
|
!
|
||||||
! that could be repeated to this determinant.
|
! that could be repeated to this determinant.
|
||||||
!
|
!
|
||||||
@ -280,16 +279,16 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
|
|||||||
enddo
|
enddo
|
||||||
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
||||||
h = h + accu_e_corr
|
h = h + accu_e_corr
|
||||||
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
|
delta_E = 1.d0/(electronic_energy(1) - h)
|
||||||
|
|
||||||
|
|
||||||
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
|
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
|
||||||
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
|
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
|
||||||
|
|
||||||
do i =2,N_st
|
do i =2,N_st
|
||||||
H_pert_diag(i) = h
|
H_pert_diag(i) = h
|
||||||
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
|
if (dabs(electronic_energy(i) - h) > 1.d-6) then
|
||||||
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
|
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
|
||||||
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
|
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
|
||||||
else
|
else
|
||||||
c_pert(i) = i_H_psi_array(i)
|
c_pert(i) = i_H_psi_array(i)
|
||||||
@ -330,11 +329,11 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments)
|
|||||||
|
|
||||||
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
|
||||||
do i =1,N_st
|
do i =1,N_st
|
||||||
if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
|
if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then
|
||||||
c_pert(i) = -1.d0
|
c_pert(i) = -1.d0
|
||||||
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
|
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
|
||||||
else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
|
else if (dabs(electronic_energy(i) - h) > 1.d-6) then
|
||||||
c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
|
c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h)
|
||||||
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
|
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
|
||||||
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
||||||
else
|
else
|
||||||
@ -350,7 +349,7 @@ end
|
|||||||
|
|
||||||
SUBST [ arguments, declarations ]
|
SUBST [ arguments, declarations ]
|
||||||
|
|
||||||
det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
|
electronic_energy,det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
|
||||||
|
|
||||||
integer, intent(in) :: Nint
|
integer, intent(in) :: Nint
|
||||||
integer, intent(in) :: ndet
|
integer, intent(in) :: ndet
|
||||||
@ -359,6 +358,7 @@ det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minili
|
|||||||
integer(bit_kind), intent(in) :: det_ref (Nint,2)
|
integer(bit_kind), intent(in) :: det_ref (Nint,2)
|
||||||
integer(bit_kind), intent(in) :: det_pert(Nint,2)
|
integer(bit_kind), intent(in) :: det_pert(Nint,2)
|
||||||
double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
|
double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
|
||||||
|
double precision , intent(in) :: electronic_energy(N_st)
|
||||||
double precision , intent(out) :: c_pert(N_st)
|
double precision , intent(out) :: c_pert(N_st)
|
||||||
double precision , intent(out) :: e_2_pert(N_st)
|
double precision , intent(out) :: e_2_pert(N_st)
|
||||||
double precision, intent(out) :: H_pert_diag(N_st)
|
double precision, intent(out) :: H_pert_diag(N_st)
|
||||||
|
3
plugins/Psiref_CAS/overwrite_with_cas.irp.f
Normal file
3
plugins/Psiref_CAS/overwrite_with_cas.irp.f
Normal file
@ -0,0 +1,3 @@
|
|||||||
|
program overwrite_w_cas
|
||||||
|
call extract_ref
|
||||||
|
end
|
@ -26,6 +26,21 @@ use bitmasks
|
|||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, psi_ref_coef_inv, (psi_det_size,n_states) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! 1/psi_ref_coef
|
||||||
|
END_DOC
|
||||||
|
integer :: i, i_state
|
||||||
|
do i_state=1,N_states
|
||||||
|
do i=1,N_det_ref
|
||||||
|
psi_ref_coef_inv(i,i_state) = 1.d0/psi_ref_coef(i,i_state)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ]
|
BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ]
|
||||||
&BEGIN_PROVIDER [ double precision, psi_ref_coef_restart, (psi_det_size,n_states) ]
|
&BEGIN_PROVIDER [ double precision, psi_ref_coef_restart, (psi_det_size,n_states) ]
|
||||||
implicit none
|
implicit none
|
||||||
|
@ -0,0 +1 @@
|
|||||||
|
Bitmask Determinants
|
24
plugins/Psiref_Utils/extract_ref.irp.f
Normal file
24
plugins/Psiref_Utils/extract_ref.irp.f
Normal file
@ -0,0 +1,24 @@
|
|||||||
|
subroutine extract_ref
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Replaces the total wave function by the normalized projection on the reference
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
integer :: i,j,k
|
||||||
|
do k=1,N_states
|
||||||
|
do j=1,N_det_ref
|
||||||
|
psi_coef(j,k) = psi_ref_coef_normalized(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do j=1,N_det_ref
|
||||||
|
do k=1,N_int
|
||||||
|
psi_det(k,1,j) = psi_ref(k,1,j)
|
||||||
|
psi_det(k,2,j) = psi_ref(k,2,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
N_det = N_det_ref
|
||||||
|
|
||||||
|
call save_wavefunction
|
||||||
|
|
||||||
|
end
|
@ -14,6 +14,47 @@ use bitmasks
|
|||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, psi_ref_coef_transp, (n_states,psi_det_size) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Transposed psi_ref_coef
|
||||||
|
END_DOC
|
||||||
|
integer :: i,j
|
||||||
|
do j=1,N_det_ref
|
||||||
|
do i=1, n_states
|
||||||
|
psi_ref_coef_transp(i,j) = psi_ref_coef(j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized, (psi_det_size,n_states) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Normalized coefficients of the reference
|
||||||
|
END_DOC
|
||||||
|
integer :: i,j,k
|
||||||
|
do k=1,N_states
|
||||||
|
do j=1,N_det_ref
|
||||||
|
psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k)
|
||||||
|
enddo
|
||||||
|
call normalize(psi_ref_coef_normalized(1,k), N_det_ref)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, psi_non_ref_coef_transp, (n_states,psi_det_size) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Transposed psi_non_ref_coef
|
||||||
|
END_DOC
|
||||||
|
integer :: i,j
|
||||||
|
do j=1,N_det_non_ref
|
||||||
|
do i=1, n_states
|
||||||
|
psi_non_ref_coef_transp(i,j) = psi_non_ref_coef(j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ]
|
BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ]
|
||||||
&BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ]
|
&BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ]
|
||||||
|
@ -184,7 +184,7 @@ def ninja_ezfio_config_rule():
|
|||||||
|
|
||||||
def get_children_of_ezfio_cfg(l_module_with_ezfio_cfg):
|
def get_children_of_ezfio_cfg(l_module_with_ezfio_cfg):
|
||||||
"""
|
"""
|
||||||
From a module list of ezfio_cfg return all the stuff create by him
|
From a module list of ezfio_cfg return all the stuff created by it
|
||||||
"""
|
"""
|
||||||
config_folder = join(QP_EZFIO, "config")
|
config_folder = join(QP_EZFIO, "config")
|
||||||
|
|
||||||
|
@ -345,7 +345,7 @@ def save_ezfio_provider(path_head, dict_code_provider):
|
|||||||
path = "{0}/ezfio_interface.irp.f".format(path_head)
|
path = "{0}/ezfio_interface.irp.f".format(path_head)
|
||||||
|
|
||||||
l_output = ["! DO NOT MODIFY BY HAND",
|
l_output = ["! DO NOT MODIFY BY HAND",
|
||||||
"! Created by $QP_ROOT/scripts/ezfio_interface.py",
|
"! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py",
|
||||||
"! from file {0}/EZFIO.cfg".format(path_head),
|
"! from file {0}/EZFIO.cfg".format(path_head),
|
||||||
"\n"]
|
"\n"]
|
||||||
|
|
||||||
|
@ -22,6 +22,7 @@ BEGIN_PROVIDER [ %(type)s, %(name)s %(size)s ]
|
|||||||
|
|
||||||
logical :: has
|
logical :: has
|
||||||
PROVIDE ezfio_filename
|
PROVIDE ezfio_filename
|
||||||
|
%(test_null_size)s
|
||||||
call ezfio_has_%(ezfio_dir)s_%(ezfio_name)s(has)
|
call ezfio_has_%(ezfio_dir)s_%(ezfio_name)s(has)
|
||||||
if (has) then
|
if (has) then
|
||||||
call ezfio_get_%(ezfio_dir)s_%(ezfio_name)s(%(name)s)
|
call ezfio_get_%(ezfio_dir)s_%(ezfio_name)s(%(name)s)
|
||||||
@ -44,6 +45,7 @@ END_PROVIDER
|
|||||||
|
|
||||||
def __repr__(self):
|
def __repr__(self):
|
||||||
self.set_write()
|
self.set_write()
|
||||||
|
self.set_test_null_size()
|
||||||
for v in self.values:
|
for v in self.values:
|
||||||
if not v:
|
if not v:
|
||||||
msg = "Error : %s is not set in EZFIO.cfg" % (v)
|
msg = "Error : %s is not set in EZFIO.cfg" % (v)
|
||||||
@ -54,8 +56,19 @@ END_PROVIDER
|
|||||||
|
|
||||||
return self.data % self.__dict__
|
return self.data % self.__dict__
|
||||||
|
|
||||||
|
def set_test_null_size(self):
|
||||||
|
if "size" not in self.__dict__:
|
||||||
|
self.__dict__["size"] = ""
|
||||||
|
if self.size != "":
|
||||||
|
self.test_null_size = "if (size(%s) == 0) return\n" % ( self.name )
|
||||||
|
else:
|
||||||
|
self.test_null_size = ""
|
||||||
|
|
||||||
def set_write(self):
|
def set_write(self):
|
||||||
self.write = ""
|
self.write = ""
|
||||||
|
if "size" in self.__dict__:
|
||||||
|
return
|
||||||
|
else:
|
||||||
if self.type in self.write_correspondance:
|
if self.type in self.write_correspondance:
|
||||||
write = self.write_correspondance[self.type]
|
write = self.write_correspondance[self.type]
|
||||||
output = self.output
|
output = self.output
|
||||||
|
@ -72,6 +72,7 @@ class H_apply(object):
|
|||||||
s["params_post"] = ""
|
s["params_post"] = ""
|
||||||
|
|
||||||
self.selection_pt2 = None
|
self.selection_pt2 = None
|
||||||
|
self.energy = "CI_electronic_energy"
|
||||||
self.perturbation = None
|
self.perturbation = None
|
||||||
self.do_double_exc = do_double_exc
|
self.do_double_exc = do_double_exc
|
||||||
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
|
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
@ -299,13 +300,13 @@ class H_apply(object):
|
|||||||
self.data["keys_work"] = """
|
self.data["keys_work"] = """
|
||||||
! if(check_double_excitation)then
|
! if(check_double_excitation)then
|
||||||
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
|
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
|
||||||
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
|
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
|
||||||
"""%(pert)
|
"""%(pert,self.energy)
|
||||||
else:
|
else:
|
||||||
self.data["keys_work"] = """
|
self.data["keys_work"] = """
|
||||||
call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
|
call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
|
||||||
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
|
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
|
||||||
"""%(pert)
|
"""%(pert,self.energy)
|
||||||
|
|
||||||
|
|
||||||
self.data["finalization"] = """
|
self.data["finalization"] = """
|
||||||
|
@ -386,7 +386,11 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,N_st,n,task_id)
|
|||||||
|
|
||||||
rc = f77_zmq_recv( zmq_socket_pull, pt2(1), 8*N_st, 0)
|
rc = f77_zmq_recv( zmq_socket_pull, pt2(1), 8*N_st, 0)
|
||||||
if (rc /= 8*N_st) then
|
if (rc /= 8*N_st) then
|
||||||
print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, pt2(1,1) , 8*N_st, 0)'
|
print *, ''
|
||||||
|
print *, ''
|
||||||
|
print *, ''
|
||||||
|
print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, pt2(1) , 8*N_st, 0)'
|
||||||
|
print *, rc
|
||||||
stop 'error'
|
stop 'error'
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
@ -10,6 +10,7 @@ subroutine $subroutine($params_main)
|
|||||||
|
|
||||||
$decls_main
|
$decls_main
|
||||||
|
|
||||||
|
integer :: i
|
||||||
integer :: i_generator
|
integer :: i_generator
|
||||||
double precision :: wall_0, wall_1
|
double precision :: wall_0, wall_1
|
||||||
integer(omp_lock_kind) :: lck
|
integer(omp_lock_kind) :: lck
|
||||||
@ -37,24 +38,22 @@ subroutine $subroutine($params_main)
|
|||||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
integer(ZMQ_PTR) :: collector_thread
|
PROVIDE nproc N_states
|
||||||
external :: $subroutine_collector
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
rc = pthread_create(collector_thread, $subroutine_collector)
|
!$OMP PRIVATE(i) &
|
||||||
|
!$OMP SHARED(zmq_socket_pair,N_states, pt2, norm_pert, H_pert_diag, n, task_id) &
|
||||||
!$OMP PARALLEL DEFAULT(private)
|
!$OMP num_threads(nproc+1)
|
||||||
!$OMP TASK PRIVATE(rc)
|
i = omp_get_thread_num()
|
||||||
rc = omp_get_thread_num()
|
if (i == 0) then
|
||||||
call $subroutine_slave_inproc(rc)
|
call $subroutine_collector()
|
||||||
!$OMP END TASK
|
integer :: n, task_id
|
||||||
!$OMP TASKWAIT
|
call pull_pt2(zmq_socket_pair, pt2, norm_pert, H_pert_diag, N_states, n, task_id)
|
||||||
|
else
|
||||||
|
call $subroutine_slave_inproc(i)
|
||||||
|
endif
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
|
||||||
integer :: n, task_id
|
|
||||||
call pull_pt2(zmq_socket_pair, pt2, norm_pert, H_pert_diag, N_st, n, task_id)
|
|
||||||
|
|
||||||
rc = pthread_join(collector_thread)
|
|
||||||
|
|
||||||
call end_zmq_pair_socket(zmq_socket_pair)
|
call end_zmq_pair_socket(zmq_socket_pair)
|
||||||
call end_parallel_job(zmq_to_qp_run_socket,'$subroutine')
|
call end_parallel_job(zmq_to_qp_run_socket,'$subroutine')
|
||||||
|
|
||||||
|
@ -315,6 +315,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
double precision, intent(inout) :: u_in(dim_in,N_st)
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
||||||
double precision, intent(out) :: energies(N_st)
|
double precision, intent(out) :: energies(N_st)
|
||||||
|
|
||||||
|
integer :: sze_8
|
||||||
integer :: iter
|
integer :: iter
|
||||||
integer :: i,j,k,l,m
|
integer :: i,j,k,l,m
|
||||||
logical :: converged
|
logical :: converged
|
||||||
@ -334,6 +335,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
double precision :: to_print(2,N_st)
|
double precision :: to_print(2,N_st)
|
||||||
double precision :: cpu, wall
|
double precision :: cpu, wall
|
||||||
|
|
||||||
|
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, Wt, y, h, lambda
|
||||||
|
|
||||||
|
|
||||||
call write_time(iunit)
|
call write_time(iunit)
|
||||||
@ -362,12 +364,15 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
enddo
|
enddo
|
||||||
write(iunit,'(A)') trim(write_buffer)
|
write(iunit,'(A)') trim(write_buffer)
|
||||||
|
|
||||||
|
integer, external :: align_double
|
||||||
|
sze_8 = align_double(sze)
|
||||||
|
|
||||||
allocate( &
|
allocate( &
|
||||||
kl_pairs(2,N_st*(N_st+1)/2), &
|
kl_pairs(2,N_st*(N_st+1)/2), &
|
||||||
W(sze,N_st,davidson_sze_max), &
|
W(sze_8,N_st,davidson_sze_max), &
|
||||||
Wt(sze), &
|
Wt(sze), &
|
||||||
U(sze,N_st,davidson_sze_max), &
|
U(sze_8,N_st,davidson_sze_max), &
|
||||||
R(sze,N_st), &
|
R(sze_8,N_st), &
|
||||||
h(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
h(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
||||||
y(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
y(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
||||||
lambda(N_st*davidson_sze_max))
|
lambda(N_st*davidson_sze_max))
|
||||||
@ -381,6 +386,8 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
! ==============
|
! ==============
|
||||||
|
|
||||||
|
|
||||||
|
if (N_st > 1) then
|
||||||
|
|
||||||
k_pairs=0
|
k_pairs=0
|
||||||
do l=1,N_st
|
do l=1,N_st
|
||||||
do k=1,l
|
do k=1,l
|
||||||
@ -393,7 +400,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
!$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, &
|
!$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, &
|
||||||
!$OMP Nint,dets_in,u_in) &
|
!$OMP Nint,dets_in,u_in) &
|
||||||
!$OMP PRIVATE(k,l,kl,i)
|
!$OMP PRIVATE(k,l,kl)
|
||||||
|
|
||||||
|
|
||||||
! Orthonormalize initial guess
|
! Orthonormalize initial guess
|
||||||
@ -415,6 +422,17 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
|
|
||||||
call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze)
|
call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
overlap(1,1) = u_dot_u(U_in(1,1),sze)
|
||||||
|
double precision :: f
|
||||||
|
f = 1.d0 / dsqrt(overlap(1,1))
|
||||||
|
do i=1,sze
|
||||||
|
U_in(i,1) = U_in(i,1) * f
|
||||||
|
enddo
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
! Davidson iterations
|
! Davidson iterations
|
||||||
! ===================
|
! ===================
|
||||||
|
|
||||||
@ -473,7 +491,10 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
! Express eigenvectors of h in the determinant basis
|
! Express eigenvectors of h in the determinant basis
|
||||||
! --------------------------------------------------
|
! --------------------------------------------------
|
||||||
|
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(k,i,l,iter2) SHARED(U,W,R,y,iter,lambda,N_st,sze)
|
||||||
do k=1,N_st
|
do k=1,N_st
|
||||||
|
!$OMP DO
|
||||||
do i=1,sze
|
do i=1,sze
|
||||||
U(i,k,iter+1) = 0.d0
|
U(i,k,iter+1) = 0.d0
|
||||||
W(i,k,iter+1) = 0.d0
|
W(i,k,iter+1) = 0.d0
|
||||||
@ -484,7 +505,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
enddo
|
enddo
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
! Compute residual vector
|
! Compute residual vector
|
||||||
! -----------------------
|
! -----------------------
|
||||||
|
@ -443,7 +443,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
|||||||
|
|
||||||
integer :: exc(0:2,2,2)
|
integer :: exc(0:2,2,2)
|
||||||
integer :: degree
|
integer :: degree
|
||||||
double precision :: get_mo_bielec_integral_schwartz
|
double precision :: get_mo_bielec_integral
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
@ -468,31 +468,31 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
|||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral_schwartz( &
|
hij = phase*get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,2), &
|
exc(2,2,2), &
|
||||||
@ -510,15 +510,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
|||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -537,15 +537,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
|||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -579,7 +579,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
|
|||||||
|
|
||||||
integer,intent(out) :: exc(0:2,2,2)
|
integer,intent(out) :: exc(0:2,2,2)
|
||||||
integer,intent(out) :: degree
|
integer,intent(out) :: degree
|
||||||
double precision :: get_mo_bielec_integral_schwartz
|
double precision :: get_mo_bielec_integral
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
@ -604,31 +604,31 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
|
|||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral_schwartz( &
|
hij = phase*get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,2), &
|
exc(2,2,2), &
|
||||||
@ -646,15 +646,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
|
|||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -673,15 +673,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
|
|||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -715,7 +715,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
|
|||||||
|
|
||||||
integer :: exc(0:2,2,2)
|
integer :: exc(0:2,2,2)
|
||||||
integer :: degree
|
integer :: degree
|
||||||
double precision :: get_mo_bielec_integral_schwartz
|
double precision :: get_mo_bielec_integral
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
@ -742,31 +742,31 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
|
|||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral_schwartz( &
|
hij = phase*get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral_schwartz( &
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral_schwartz( &
|
get_mo_bielec_integral( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,2), &
|
exc(2,2,2), &
|
||||||
@ -784,15 +784,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
|
|||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -811,15 +811,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
|
|||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
if (.not.has_mipi(i)) then
|
||||||
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
@ -1608,12 +1608,11 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
|||||||
integer :: i,j,k,l, jj,ii
|
integer :: i,j,k,l, jj,ii
|
||||||
integer :: i0, j0
|
integer :: i0, j0
|
||||||
|
|
||||||
integer, allocatable :: shortcut(:), sort_idx(:)
|
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
|
||||||
integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
|
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
|
||||||
integer(bit_kind) :: sorted_i(Nint)
|
integer(bit_kind) :: sorted_i(Nint)
|
||||||
|
|
||||||
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi
|
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi
|
||||||
double precision :: local_threshold
|
|
||||||
|
|
||||||
|
|
||||||
ASSERT (Nint > 0)
|
ASSERT (Nint > 0)
|
||||||
@ -1621,104 +1620,83 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
|||||||
ASSERT (n>0)
|
ASSERT (n>0)
|
||||||
PROVIDE ref_bitmask_energy davidson_criterion
|
PROVIDE ref_bitmask_energy davidson_criterion
|
||||||
|
|
||||||
allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n))
|
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
|
||||||
v_0 = 0.d0
|
v_0 = 0.d0
|
||||||
|
|
||||||
call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
|
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
|
||||||
|
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)&
|
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i)&
|
||||||
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version)
|
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version)
|
||||||
allocate(vt(n))
|
allocate(vt(n))
|
||||||
Vt = 0.d0
|
Vt = 0.d0
|
||||||
|
|
||||||
!$OMP DO SCHEDULE(dynamic)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
do sh=1,shortcut(0,1)
|
||||||
do sh2=1,sh
|
do sh2=sh,shortcut(0,1)
|
||||||
exa = 0
|
exa = 0
|
||||||
do ni=1,Nint
|
do ni=1,Nint
|
||||||
exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2)))
|
exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
|
||||||
end do
|
end do
|
||||||
if(exa > 2) then
|
if(exa > 2) then
|
||||||
cycle
|
cycle
|
||||||
end if
|
end if
|
||||||
|
|
||||||
do i=shortcut(sh),shortcut(sh+1)-1
|
do i=shortcut(sh,1),shortcut(sh+1,1)-1
|
||||||
org_i = sort_idx(i)
|
org_i = sort_idx(i,1)
|
||||||
local_threshold = threshold_davidson - dabs(u_0(org_i))
|
|
||||||
if(sh==sh2) then
|
if(sh==sh2) then
|
||||||
endi = i-1
|
endi = i-1
|
||||||
else
|
else
|
||||||
endi = shortcut(sh2+1)-1
|
endi = shortcut(sh2+1,1)-1
|
||||||
end if
|
end if
|
||||||
do ni=1,Nint
|
do ni=1,Nint
|
||||||
sorted_i(ni) = sorted(ni,i)
|
sorted_i(ni) = sorted(ni,i,1)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do j=shortcut(sh2),endi
|
do j=shortcut(sh2,1),endi
|
||||||
org_j = sort_idx(j)
|
org_j = sort_idx(j,1)
|
||||||
if ( dabs(u_0(org_j)) > local_threshold ) then
|
|
||||||
ext = exa
|
ext = exa
|
||||||
do ni=1,Nint
|
do ni=1,Nint
|
||||||
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j)))
|
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
|
||||||
end do
|
end do
|
||||||
if(ext <= 4) then
|
if(ext <= 4) then
|
||||||
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
|
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
|
||||||
vt (org_i) = vt (org_i) + hij*u_0(org_j)
|
vt (org_i) = vt (org_i) + hij*u_0(org_j)
|
||||||
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
||||||
endif
|
endif
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
!$OMP CRITICAL
|
|
||||||
do i=1,n
|
|
||||||
v_0(i) = v_0(i) + vt(i)
|
|
||||||
enddo
|
|
||||||
!$OMP END CRITICAL
|
|
||||||
|
|
||||||
deallocate(vt)
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
|
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)&
|
|
||||||
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version)
|
|
||||||
allocate(vt(n))
|
|
||||||
Vt = 0.d0
|
|
||||||
|
|
||||||
!$OMP DO SCHEDULE(dynamic)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
do sh=1,shortcut(0,2)
|
||||||
do i=shortcut(sh),shortcut(sh+1)-1
|
do i=shortcut(sh,2),shortcut(sh+1,2)-1
|
||||||
org_i = sort_idx(i)
|
org_i = sort_idx(i,2)
|
||||||
local_threshold = threshold_davidson - dabs(u_0(org_i))
|
do j=shortcut(sh,2),i-1
|
||||||
do j=shortcut(sh),i-1
|
org_j = sort_idx(j,2)
|
||||||
org_j = sort_idx(j)
|
|
||||||
if ( dabs(u_0(org_j)) > local_threshold ) then
|
|
||||||
ext = 0
|
ext = 0
|
||||||
do ni=1,Nint
|
do ni=1,Nint
|
||||||
ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j)))
|
ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
|
||||||
end do
|
end do
|
||||||
if(ext == 4) then
|
if(ext == 4) then
|
||||||
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
|
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
|
||||||
vt (org_i) = vt (org_i) + hij*u_0(org_j)
|
vt (org_i) = vt (org_i) + hij*u_0(org_j)
|
||||||
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
||||||
end if
|
end if
|
||||||
end if
|
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
!$OMP CRITICAL
|
!$OMP CRITICAL
|
||||||
do i=1,n
|
do i=n,1,-1
|
||||||
v_0(i) = v_0(i) + vt(i)
|
v_0(i) = v_0(i) + vt(i)
|
||||||
enddo
|
enddo
|
||||||
!$OMP END CRITICAL
|
!$OMP END CRITICAL
|
||||||
|
|
||||||
deallocate(vt)
|
deallocate(vt)
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
@ -1728,3 +1706,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
|||||||
deallocate (shortcut, sort_idx, sorted, version)
|
deallocate (shortcut, sort_idx, sorted, version)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine apply_excitation(det, exc, res, ok, Nint)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: Nint
|
||||||
|
integer, intent(in) :: exc(0:2,2,2)
|
||||||
|
integer(bit_kind),intent(in) :: det(Nint, 2)
|
||||||
|
integer(bit_kind),intent(out) :: res(Nint, 2)
|
||||||
|
logical, intent(out) :: ok
|
||||||
|
integer :: h1,p1,h2,p2,s1,s2,degree
|
||||||
|
integer :: ii, pos
|
||||||
|
|
||||||
|
|
||||||
|
ok = .false.
|
||||||
|
degree = exc(0,1,1) + exc(0,1,2)
|
||||||
|
|
||||||
|
if(.not. (degree > 0 .and. degree <= 2)) then
|
||||||
|
print *, degree
|
||||||
|
print *, "apply ex"
|
||||||
|
STOP
|
||||||
|
endif
|
||||||
|
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
|
res = det
|
||||||
|
|
||||||
|
ii = (h1-1)/bit_kind_size + 1
|
||||||
|
pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
|
||||||
|
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
|
||||||
|
res(ii, s1) = ibclr(res(ii, s1), pos)
|
||||||
|
|
||||||
|
ii = (p1-1)/bit_kind_size + 1
|
||||||
|
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
|
||||||
|
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
|
||||||
|
res(ii, s1) = ibset(res(ii, s1), pos)
|
||||||
|
|
||||||
|
if(degree == 2) then
|
||||||
|
ii = (h2-1)/bit_kind_size + 1
|
||||||
|
pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
|
||||||
|
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
|
||||||
|
res(ii, s2) = ibclr(res(ii, s2), pos)
|
||||||
|
|
||||||
|
ii = (p2-1)/bit_kind_size + 1
|
||||||
|
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
|
||||||
|
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
|
||||||
|
res(ii, s2) = ibset(res(ii, s2), pos)
|
||||||
|
endif
|
||||||
|
|
||||||
|
ok = .true.
|
||||||
|
end subroutine
|
||||||
|
|
||||||
|
|
||||||
|
@ -374,20 +374,16 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
|
|||||||
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
integer(ZMQ_PTR) :: collector_thread
|
PROVIDE nproc
|
||||||
external :: ao_bielec_integrals_in_map_collector
|
!$OMP PARALLEL DEFAULT(private) num_threads(nproc+1)
|
||||||
rc = pthread_create(collector_thread, ao_bielec_integrals_in_map_collector)
|
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(private)
|
|
||||||
!$OMP TASK PRIVATE(i)
|
|
||||||
i = omp_get_thread_num()
|
i = omp_get_thread_num()
|
||||||
|
if (i==0) then
|
||||||
|
call ao_bielec_integrals_in_map_collector(i)
|
||||||
|
else
|
||||||
call ao_bielec_integrals_in_map_slave_inproc(i)
|
call ao_bielec_integrals_in_map_slave_inproc(i)
|
||||||
!$OMP END TASK
|
endif
|
||||||
!$OMP TASKWAIT
|
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
rc = pthread_join(collector_thread)
|
|
||||||
|
|
||||||
call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals')
|
call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals')
|
||||||
|
|
||||||
|
|
||||||
|
@ -67,6 +67,8 @@ end
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
|
subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
|
||||||
use map_module
|
use map_module
|
||||||
use f77_zmq
|
use f77_zmq
|
||||||
@ -107,7 +109,7 @@ subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
|
|||||||
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, 0)
|
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, 0)
|
||||||
enddo
|
enddo
|
||||||
call compute_ao_integrals_jl(l,l,n_integrals,buffer_i,buffer_value)
|
call compute_ao_integrals_jl(l,l,n_integrals,buffer_i,buffer_value)
|
||||||
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_integrals)
|
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
|
||||||
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
|
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
@ -127,7 +129,7 @@ subroutine pull_integrals(zmq_socket_pull, n_integrals, buffer_i, buffer_value,
|
|||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! How the collector pulls the computed integrals
|
! How the collector pulls the computed integrals
|
||||||
END_DOC
|
END_DOC
|
||||||
integer(ZMQ_PTR), intent(out) :: zmq_socket_pull
|
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||||
integer, intent(out) :: n_integrals
|
integer, intent(out) :: n_integrals
|
||||||
integer(key_kind), intent(out) :: buffer_i(*)
|
integer(key_kind), intent(out) :: buffer_i(*)
|
||||||
real(integral_kind), intent(out) :: buffer_value(*)
|
real(integral_kind), intent(out) :: buffer_value(*)
|
||||||
|
@ -323,9 +323,9 @@ double precision function mo_bielec_integral(i,j,k,l)
|
|||||||
! Returns one integral <ij|kl> in the MO basis
|
! Returns one integral <ij|kl> in the MO basis
|
||||||
END_DOC
|
END_DOC
|
||||||
integer, intent(in) :: i,j,k,l
|
integer, intent(in) :: i,j,k,l
|
||||||
double precision :: get_mo_bielec_integral_schwartz
|
double precision :: get_mo_bielec_integral
|
||||||
PROVIDE mo_bielec_integrals_in_map
|
PROVIDE mo_bielec_integrals_in_map
|
||||||
mo_bielec_integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map)
|
mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -8,7 +8,8 @@ program qp_ao_ints
|
|||||||
|
|
||||||
call switch_qp_run_to_master
|
call switch_qp_run_to_master
|
||||||
|
|
||||||
PROVIDE zmq_context
|
zmq_context = f77_zmq_ctx_new ()
|
||||||
|
|
||||||
! Set the state of the ZMQ
|
! Set the state of the ZMQ
|
||||||
zmq_state = 'ao_integrals'
|
zmq_state = 'ao_integrals'
|
||||||
|
|
||||||
|
@ -3,10 +3,14 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
|
|||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Pseudo-potential
|
! Pseudo-potential
|
||||||
END_DOC
|
END_DOC
|
||||||
if (do_pseudo) then
|
|
||||||
ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local
|
|
||||||
else
|
|
||||||
ao_pseudo_integral = 0.d0
|
ao_pseudo_integral = 0.d0
|
||||||
|
if (do_pseudo) then
|
||||||
|
if (pseudo_klocmax > 0) then
|
||||||
|
ao_pseudo_integral += ao_pseudo_integral_local
|
||||||
|
endif
|
||||||
|
if (pseudo_kmax > 0) then
|
||||||
|
ao_pseudo_integral += ao_pseudo_integral_non_local
|
||||||
|
endif
|
||||||
endif
|
endif
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
14
src/MO_Basis/swap_mos.irp.f
Normal file
14
src/MO_Basis/swap_mos.irp.f
Normal file
@ -0,0 +1,14 @@
|
|||||||
|
program swap_mos
|
||||||
|
implicit none
|
||||||
|
integer :: i,j, i1, i2
|
||||||
|
double precision :: x
|
||||||
|
print *, 'MOs to swap?'
|
||||||
|
read(*,*) i1, i2
|
||||||
|
do i=1,ao_num_align
|
||||||
|
x = mo_coef(i,i1)
|
||||||
|
mo_coef(i,i1) = mo_coef(i,i2)
|
||||||
|
mo_coef(i,i2) = x
|
||||||
|
enddo
|
||||||
|
call save_mos
|
||||||
|
|
||||||
|
end
|
@ -73,6 +73,10 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
|
|||||||
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
|
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
|
||||||
integer :: info, i, j
|
integer :: info, i, j
|
||||||
|
|
||||||
|
if (n < 2) then
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n))
|
allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n))
|
||||||
|
|
||||||
call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
|
call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
|
||||||
@ -151,6 +155,10 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
|
|||||||
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
|
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
|
||||||
integer :: info, i, j, k
|
integer :: info, i, j, k
|
||||||
|
|
||||||
|
if (n < 2) then
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
call svd(overlap,lda,U,ldc,D,Vt,lda,m,n)
|
call svd(overlap,lda,U,ldc,D,Vt,lda,m,n)
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
@ -295,18 +295,6 @@ BEGIN_PROVIDER [ integer, nproc ]
|
|||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ integer, iproc_save, (nproc) ]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! iproc_save(i) = i-1. Used to start threads with pthreads.
|
|
||||||
END_DOC
|
|
||||||
integer :: i
|
|
||||||
do i=1,nproc
|
|
||||||
iproc_save(i) = i-1
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
double precision function u_dot_v(u,v,sze)
|
double precision function u_dot_v(u,v,sze)
|
||||||
implicit none
|
implicit none
|
||||||
@ -324,6 +312,7 @@ double precision function u_dot_v(u,v,sze)
|
|||||||
t3 = t2+t2
|
t3 = t2+t2
|
||||||
t4 = t3+t2
|
t4 = t3+t2
|
||||||
u_dot_v = 0.d0
|
u_dot_v = 0.d0
|
||||||
|
!DIR$ VECTOR ALWAYS
|
||||||
do i=1,t2
|
do i=1,t2
|
||||||
u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) + &
|
u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) + &
|
||||||
u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i)
|
u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i)
|
||||||
@ -359,6 +348,7 @@ double precision function u_dot_u(u,sze)
|
|||||||
! u_dot_u = u_dot_u+u(i)*u(i)
|
! u_dot_u = u_dot_u+u(i)*u(i)
|
||||||
! enddo
|
! enddo
|
||||||
|
|
||||||
|
!DIR$ VECTOR ALWAYS
|
||||||
do i=1,sze
|
do i=1,sze
|
||||||
u_dot_u = u_dot_u + u(i)*u(i)
|
u_dot_u = u_dot_u + u(i)*u(i)
|
||||||
enddo
|
enddo
|
||||||
|
@ -361,6 +361,8 @@ subroutine end_zmq_pull_socket(zmq_socket_pull)
|
|||||||
stop 'error'
|
stop 'error'
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
call sleep(1) ! see https://github.com/zeromq/libzmq/issues/1922
|
||||||
|
|
||||||
rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4)
|
rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4)
|
||||||
if (rc /= 0) then
|
if (rc /= 0) then
|
||||||
stop 'Unable to set ZMQ_LINGER on zmq_socket_pull'
|
stop 'Unable to set ZMQ_LINGER on zmq_socket_pull'
|
||||||
|
@ -53,10 +53,10 @@ function test_exe() {
|
|||||||
}
|
}
|
||||||
|
|
||||||
function run_HF() {
|
function run_HF() {
|
||||||
thresh=1.e-8
|
thresh=1.e-7
|
||||||
test_exe SCF || skip
|
test_exe SCF || skip
|
||||||
ezfio set_file $1
|
ezfio set_file $1
|
||||||
ezfio set hartree_fock thresh_scf 1.e-10
|
ezfio set hartree_fock thresh_scf 1.e-11
|
||||||
qp_run SCF $1
|
qp_run SCF $1
|
||||||
energy="$(ezfio get hartree_fock energy)"
|
energy="$(ezfio get hartree_fock energy)"
|
||||||
eq $energy $2 $thresh
|
eq $energy $2 $thresh
|
||||||
@ -155,7 +155,7 @@ function run_all_1h_1p() {
|
|||||||
ezfio set determinants read_wf True
|
ezfio set determinants read_wf True
|
||||||
qp_run mrcc_cassd $INPUT
|
qp_run mrcc_cassd $INPUT
|
||||||
energy="$(ezfio get mrcc_cassd energy)"
|
energy="$(ezfio get mrcc_cassd energy)"
|
||||||
eq $energy -0.762303253805911E+02 1.E-3
|
eq $energy -76.2288648023833 1.e-4
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -166,7 +166,7 @@ function run_all_1h_1p() {
|
|||||||
}
|
}
|
||||||
|
|
||||||
@test "SCF H2O VDZ pseudo" {
|
@test "SCF H2O VDZ pseudo" {
|
||||||
run_HF h2o_pseudo.ezfio -0.169483703904991E+02
|
run_HF h2o_pseudo.ezfio -16.9483703905461
|
||||||
}
|
}
|
||||||
|
|
||||||
@test "FCI H2O VDZ pseudo" {
|
@test "FCI H2O VDZ pseudo" {
|
||||||
|
Loading…
Reference in New Issue
Block a user