LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-25 05:43:47 +01:00
Code Issues Releases Wiki Activity

merged with scemama

Browse Source
This commit is contained in:
Yann Garniron 2018-02-09 09:58:19 +01:00
commit ff1ef59d32
67 changed files with 643 additions and 2041 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.travis.yml
View File

@ -24,6 +24,7 @@ addons:
cache:
directories:
- $HOME/.opam/
- $HOME/lapack-release
language: python
python:

6
README.md
View File

@ -1,7 +1,9 @@
## IMPORTANT
If you have problems upgrading to the current version, consider re-installing everything from scratch including the OCaml compiler.
To do this, you will have to remove the `quantum_package` directory **and** the `$HOME/.opam` directory as well.
If you have problems upgrading to the current version, first try
`qp_upgrade_ocaml.sh`. If it fails, then consider re-installing everything from
scratch including the OCaml compiler. To do this, you will have to remove the
`quantum_package` directory **and** the `$HOME/.opam` directory as well.

28
configure vendored
View File

@ -49,7 +49,7 @@ QP_ROOT_INSTALL = join(QP_ROOT, "install")
os.environ["PATH"] = os.environ["PATH"] + ":" + QP_ROOT_BIN
d_dependency = {
"ocaml": ["m4", "curl", "zlib", "patch", "gcc", "zeromq"],
"ocaml": ["m4", "curl", "zlib", "patch", "gcc", "zeromq", "gmp"],
"m4": ["make"],
"curl": ["make"],
"zlib": ["gcc", "make"],
@ -67,7 +67,8 @@ d_dependency = {
"ninja": ["g++", "python"],
"make": [],
"p_graphviz": ["python"],
"bats": []
"bats": [],
"gmp" : ["make", "g++"]
}
from collections import namedtuple
@ -136,6 +137,11 @@ zeromq = Info(
description=' ZeroMQ',
default_path=join(QP_ROOT_LIB, "libzmq.a"))
gmp= Info(
url='https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2',
description=' The GNU Multiple Precision Arithmetic Library',
default_path=join(QP_ROOT_LIB, "libgmp.a"))
f77zmq = Info(
url='{head}/zeromq/f77_zmq/{tail}'.format(**path_github),
description=' F77-ZeroMQ',
@ -155,7 +161,7 @@ d_info = dict()
for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt",
"resultsFile", "ninja", "emsl", "ezfio", "p_graphviz",
"zeromq", "f77zmq", "bats"]:
"zeromq", "f77zmq", "bats", "gmp"]:
exec ("d_info['{0}']={0}".format(m))
@ -480,16 +486,16 @@ def create_ninja_and_rc(l_installed):
'export QP_PYTHON={0}'.format(":".join(l_python)), "",
'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
'function qp_append_export () {',
' #Append path $2:${!1}. Add the semicolon only if ${!1} is defined',
'function qp_prepend_export () {',
' #Prepend path $2:${!1}. Add the semicolon only if ${!1} is defined',
' eval "value_1=\"\${$1}\""',
' echo ${2}${value_1:+:${value_1}}',
' echo ${value_1:+${value_1}:}${2}',
'}',
'export PYTHONPATH=$(qp_append_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")',
'export PATH=$(qp_append_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)',
'export LD_LIBRARY_PATH=$(qp_append_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
'export LIBRARY_PATH=$(qp_append_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
'export C_INCLUDE_PATH=$(qp_append_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)',
'export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")',
'export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)',
'export LD_LIBRARY_PATH=$(qp_prepend_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
'export LIBRARY_PATH=$(qp_prepend_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
'export C_INCLUDE_PATH=$(qp_prepend_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)',
'',
'if [[ $SHELL == "bash" ]] ; then',
' source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh',

2
install/scripts/install_lapack.sh
View File

@ -1,6 +1,6 @@
#!/bin/bash -x
git clone https://github.com/Reference-LAPACK/lapack-release.git
git clone https://github.com/Reference-LAPACK/lapack-release.git || echo "Clone failed"
cd lapack-release
cp make.inc.example make.inc
make -j 8

5
install/scripts/install_ocaml.sh
View File

@ -5,11 +5,12 @@ QP_ROOT=$PWD
cd -
# Normal installation
PACKAGES="core.v0.9.1 cryptokit.1.10 ocamlfind sexplib.v0.9.1 ZMQ ppx_sexp_conv ppx_deriving"
PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 ZMQ ppx_sexp_conv ppx_deriving"
# Needed for ZeroMQ
export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}"
export LIBRARY_PATH="${QP_ROOT}"/lib:"${LIBRARY_PATH}"
export LDFLAGS="-L$QP_ROOT/lib"
export LD_LIBRARY_PATH="${QP_ROOT}"/lib:"${LD_LIBRARY_PATH}"
# return 0 if program version is equal or greater than check version
@ -64,7 +65,7 @@ fi
cd Downloads || exit 1
chmod +x ocaml.sh || exit 1
echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.04.2 || exit 1
echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.06.0 || exit 1
${QP_ROOT}/bin/opam config setup -a -q || exit 1

12
ocaml/Input_determinants_by_hand.ml
View File

@ -93,8 +93,16 @@ end = struct
;;
let write_n_states n =
States_number.to_int n
|> Ezfio.set_determinants_n_states
let n_states =
States_number.to_int n
in
Ezfio.set_determinants_n_states n_states;
let data =
Array.create n_states 1.
|> Array.to_list
in
Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_states |] ~data
|> Ezfio.set_determinants_state_average_weight
;;
let write_state_average_weight data =

11
ocaml/TaskServer.ml
View File

@ -1,6 +1,7 @@
open Core
open Qptypes
module StringHashtbl = Hashtbl.Make(String)
type pub_state =
| Waiting
@ -28,7 +29,7 @@ type t =
progress_bar : Progress_bar.t option ;
running : bool;
accepting_clients : bool;
data : (string, string) Hashtbl.t;
data : string StringHashtbl.t;
}
@ -208,7 +209,7 @@ let end_job msg program_state rep_socket pair_socket =
address_inproc = None;
running = true;
accepting_clients = false;
data = Hashtbl.create ~hashable:String.hashable ();
data = StringHashtbl.create ();
}
and wait n =
@ -592,7 +593,7 @@ let put_data msg rest_of_msg program_state rep_socket =
in
let success () =
Hashtbl.set program_state.data ~key ~data:value ;
StringHashtbl.set program_state.data ~key ~data:value ;
Message.PutDataReply (Message.PutDataReply_msg.create ())
|> Message.to_string
|> ZMQ.Socket.send rep_socket;
@ -622,7 +623,7 @@ let get_data msg program_state rep_socket =
let success () =
let value =
match Hashtbl.find program_state.data key with
match StringHashtbl.find program_state.data key with
| Some value -> value
| None -> ""
in
@ -776,7 +777,7 @@ let run ~port =
address_inproc = None;
progress_bar = None ;
accepting_clients = false;
data = Hashtbl.create ~hashable:String.hashable ();
data = StringHashtbl.create ();
}
in

2
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -665,7 +665,7 @@ let run ?o b au c d m p cart xyz_file =
let command =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:(fun () -> "

2
ocaml/qp_create_guess.ml
View File

@ -128,7 +128,7 @@ let spec =
+> anon ("ezfio_file" %: string)
let () =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:( fun () -> "
Creates an open-shell multiplet initial guess\n\n" )

2
ocaml/qp_find_pi_space.ml
View File

@ -95,7 +95,7 @@ let spec =
let command =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:(fun () ->
"Find all the pi molecular orbitals to create a pi space.

2
ocaml/qp_print.ml
View File

@ -141,7 +141,7 @@ let run_o ~action ezfio_filename =
;;
let command =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:(fun () ->
"

2
ocaml/qp_run.ml
View File

@ -150,7 +150,7 @@ let spec =
let () =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:( fun () -> "
Executes a Quantum Package binary file among these:\n\n"

2
ocaml/qp_set_mo_class.ml
View File

@ -323,7 +323,7 @@ let spec =
let command =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:(fun () ->
"Set the orbital classes in an EZFIO directory

2
plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES
View File

@ -1,2 +1,2 @@
Generators_CAS Perturbation Selectors_CASSD ZMQ
Generators_CAS Perturbation Selectors_CASSD ZMQ DavidsonUndressed

2
plugins/CID/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod Davidson
Selectors_full SingleRefMethod DavidsonUndressed

2
plugins/CID_selected/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation CID
Perturbation CID DavidsonUndressed

2
plugins/CIS/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod Davidson
Selectors_full SingleRefMethod DavidsonUndressed

2
plugins/CISD/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod Davidson
Selectors_full SingleRefMethod DavidsonUndressed

2
plugins/Casino/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants Davidson
Determinants DavidsonUndressed

2
plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_full ZMQ FourIdx MPI
Perturbation Selectors_full Generators_full ZMQ FourIdx MPI DavidsonUndressed

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

@ -25,8 +25,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
double precision, external :: omp_get_wtime
double precision :: state_average_weight_save(N_states), w(N_states)
double precision :: time
double precision :: w(N_states)
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
if (N_det < max(10,N_states)) then
@ -35,18 +35,19 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
error(:) = 0.d0
else
state_average_weight_save(:) = state_average_weight(:)
do pt2_stoch_istate=1,N_states
SOFT_TOUCH pt2_stoch_istate
w(:) = 0.d0
w(pt2_stoch_istate) = 1.d0
call update_psi_average_norm_contrib(w)
state_average_weight(:) = 0.d0
state_average_weight(pt2_stoch_istate) = 1.d0
TOUCH state_average_weight
allocate(pt2_detail(N_states,N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
sumabove = 0d0
sum2above = 0d0
Nabove = 0d0
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors
computed = .false.
@ -141,7 +142,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
deallocate(pt2_detail, comb, computed, tbc)
enddo
FREE psi_average_norm_contrib pt2_stoch_istate
FREE pt2_stoch_istate
state_average_weight(:) = state_average_weight_save(:)
TOUCH state_average_weight
endif
do k=N_det+1,N_states
pt2(k) = 0.d0

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

@ -623,6 +623,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
min_e_pert = 0d0
! double precision :: hij
! call i_h_j(psi_det_generators(1,1,i_generator), det, N_int, hij)
do istate=1,N_states
delta_E = E0(istate) - Hii
val = mat(istate, p1, p2) + mat(istate, p1, p2)
@ -633,7 +636,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
e_pert = 0.5d0 * (tmp - delta_E)
pt2(istate) = pt2(istate) + e_pert
min_e_pert = min(e_pert,min_e_pert)
! ci(istate) = e_pert / mat(istate, p1, p2)
! ci(istate) = e_pert / hij
end do
if(min_e_pert <= buf%mini) then

1
plugins/GPI2/NEEDED_CHILDREN_MODULES
View File

@ -1 +0,0 @@
Determinants

13
plugins/GPI2/gpi_test.irp.f
View File

@ -1,13 +0,0 @@
program test
double precision :: energy(N_states)
if (is_gaspi_master) then
energy = 1.d0
else
energy = 0.d0
endif
call broadcast_wf(energy)
print *, 'energy (1.d0) :', GASPI_rank, energy(1)
print *, 'coef :', GASPI_rank, psi_coef(1,1)
print *, 'det :', GASPI_rank, psi_det (1,1,1)
call gaspi_finalize
end

76
plugins/GPI2/utils.irp.f
View File

@ -1,76 +0,0 @@
BEGIN_PROVIDER [ logical, GASPI_is_initialized ]
&BEGIN_PROVIDER [ logical, has_gaspi ]
implicit none
BEGIN_DOC
! This is true when GASPI_Init has been called
END_DOC
has_gaspi = .False.
IRP_IF GASPI
use GASPI
integer(gaspi_return_t) :: res
res = gaspi_proc_init(GASPI_BLOCK)
if (res /= GASPI_SUCCESS) then
print *, res
print *, 'GASPI failed to initialize'
stop -1
endif
has_gaspi = .True.
IRP_ENDIF
GASPI_is_initialized = .True.
END_PROVIDER
BEGIN_PROVIDER [ integer, GASPI_rank ]
&BEGIN_PROVIDER [ integer, GASPI_size ]
&BEGIN_PROVIDER [ logical, is_GASPI_master ]
implicit none
BEGIN_DOC
! Usual GASPI variables
END_DOC
PROVIDE GASPI_is_initialized
IRP_IF GASPI
use GASPI
integer(gaspi_return_t) :: res
integer(gaspi_rank_t) :: n
res = gaspi_proc_num(n)
GASPI_size = n
if (res /= GASPI_SUCCESS) then
print *, res
print *, 'Unable to get GASPI_size'
stop -1
endif
res = gaspi_proc_rank(n)
GASPI_rank = n
if (res /= GASPI_SUCCESS) then
print *, res
print *, 'Unable to get GASPI_rank'
stop -1
endif
is_GASPI_master = (GASPI_rank == 0)
IRP_ELSE
GASPI_rank = 0
GASPI_size = 1
is_GASPI_master = .True.
IRP_ENDIF
END_PROVIDER
subroutine gaspi_finalize()
implicit none
PROVIDE GASPI_is_initialized
IRP_IF GASPI
use GASPI
integer(gaspi_return_t) :: res
res = gaspi_proc_term(GASPI_BLOCK)
if (res /= GASPI_SUCCESS) then
print *, res
print *, 'Unable to finalize GASPI'
stop -1
endif
IRP_ENDIF
end subroutine

11
plugins/Generators_full/generators.irp.f
View File

@ -30,15 +30,8 @@ END_PROVIDER
! Hartree-Fock determinant
END_DOC
integer :: i, k
psi_coef_generators = 0.d0
psi_det_generators = 0_bit_kind
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_det_sorted(k,1,i)
psi_det_generators(k,2,i) = psi_det_sorted(k,2,i)
enddo
psi_coef_generators(i,:) = psi_coef_sorted(i,:)
enddo
psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
END_PROVIDER

2
plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
View File

@ -14,6 +14,8 @@ END_DOC
integer :: i,j
double precision, allocatable :: mo_coef_save(:,:)
PROVIDE ao_md5 mo_occ level_shift
allocate(mo_coef_save(ao_num,mo_tot_num), &
Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS), &

2
plugins/Hartree_Fock/SCF.irp.f
View File

@ -23,7 +23,7 @@ subroutine create_guess
mo_coef = ao_ortho_lowdin_coef
TOUCH mo_coef
mo_label = 'Guess'
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,1,.false.)
SOFT_TOUCH mo_coef mo_label
else if (mo_guess_type == "Huckel") then
call huckel_guess

10
plugins/MRCC_Utils/H_apply.irp.f
View File

@ -3,19 +3,17 @@ BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("mrcc")
s.data["parameters"] = ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["parameters"] = ", delta_ij_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["declarations"] += """
integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, Ndet_ref)
double precision, intent(in) :: delta_ii_(Nstates, Ndet_ref)
"""
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_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["keys_work"] = "call mrcc_dress(delta_ij_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
s.data["params_post"] += ", delta_ij_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["params_main"] += "delta_ij_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["decls_main"] += """
integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates
double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref)
"""
s.data["finalization"] = ""
s.data["copy_buffer"] = ""

1128
plugins/MRCC_Utils/davidson.irp.f
View File

File diff suppressed because it is too large Load Diff

5
plugins/MRCC_Utils/mrcc_dress.irp.f
View File

@ -14,14 +14,13 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
END_PROVIDER
subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
subroutine mrcc_dress(delta_ij_, Nstates, Ndet_non_ref, Ndet_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) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l,m
@ -265,10 +264,8 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
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)
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)
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) + 0.5d0 * dIa_hla(i_state,k_sd)

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

@ -139,210 +139,6 @@ BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
implicit none
BEGIN_DOC
! Dressed H with Delta_ij
END_DOC
integer :: i, j,istate,ii,jj
do istate = 1,N_states
do j=1,N_det
do i=1,N_det
h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j)
enddo
enddo
do ii = 1, N_det_ref
i =idx_ref(ii)
h_matrix_dressed(i,i,istate) += delta_ii(istate,ii)
do jj = 1, N_det_non_ref
j =idx_non_ref(jj)
h_matrix_dressed(i,j,istate) += delta_ij(istate,jj,ii)
h_matrix_dressed(j,i,istate) += delta_ij(istate,jj,ii)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the dressed CI matrix
END_DOC
double precision :: ovrlp,u_dot_v
integer :: i_good_state
integer, allocatable :: index_good_state_array(:)
logical, allocatable :: good_state_array(:)
double precision, allocatable :: s2_values_tmp(:)
integer :: i_other_state
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
integer :: i_state
double precision :: e_0
integer :: i,j,k
double precision, allocatable :: s2_eigvalues(:)
double precision, allocatable :: e_array(:)
integer, allocatable :: iorder(:)
integer :: mrcc_state
do j=1,min(N_states,N_det)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = psi_coef(i,j)
enddo
enddo
if (diag_algorithm == "Davidson") then
allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),&
eigenvalues(size(CI_electronic_energy_dressed,1)))
do j=1,min(N_states,N_det)
do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
do mrcc_state=1,N_states
do j=mrcc_state,min(N_states,N_det)
do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
call davidson_diag_mrcc_HS2(psi_det,eigenvectors, &
size(eigenvectors,1), &
eigenvalues,N_det,N_states,N_states_diag,N_int, &
6,mrcc_state)
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo
do k=N_states+1,N_states_diag
CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
CI_electronic_energy_dressed(k) = eigenvalues(k)
enddo
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors_dressed,1))
deallocate (eigenvectors,eigenvalues)
else if (diag_algorithm == "Lapack") then
allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_dressed,size(H_matrix_dressed,1),N_det)
CI_electronic_energy_dressed(:) = 0.d0
if (s2_eig) then
i_state = 0
allocate (s2_eigvalues(N_det))
allocate(index_good_state_array(N_det),good_state_array(N_det))
good_state_array = .False.
call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
N_det,size(eigenvectors,1))
do j=1,N_det
! Select at least n_states states with S^2 values closed to "expected_s2"
if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
i_state += 1
index_good_state_array(i_state) = j
good_state_array(j) = .True.
endif
if (i_state==N_states) then
exit
endif
enddo
if (i_state /= 0) then
! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j))
enddo
CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j))
CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j))
enddo
i_other_state = 0
do j = 1, N_det
if(good_state_array(j))cycle
i_other_state +=1
if(i_state+i_other_state.gt.n_states_diag)then
exit
endif
do i=1,N_det
CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j)
CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
enddo
else
print*,''
print*,'!!!!!!!! WARNING !!!!!!!!!'
print*,' Within the ',N_det,'determinants selected'
print*,' and the ',N_states_diag,'states requested'
print*,' We did not find any state with S^2 values close to ',expected_s2
print*,' We will then set the first N_states eigenvectors of the H matrix'
print*,' as the CI_eigenvectors_dressed'
print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
print*,''
do j=1,min(N_states_diag,N_det)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(j) = eigenvalues(j)
CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j)
enddo
endif
deallocate(index_good_state_array,good_state_array)
deallocate(s2_eigvalues)
else
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states_diag),size(eigenvectors,1))
! Select the "N_states_diag" states of lowest energy
do j=1,min(N_det,N_states_diag)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(j) = eigenvalues(j)
enddo
endif
deallocate(eigenvectors,eigenvalues)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the dressed CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(6)
do j=1,min(N_det,N_states)
write(st,'(I4)') j
CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st))
call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
enddo
END_PROVIDER
subroutine diagonalize_CI_dressed(lambda)
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
double precision, intent(in) :: lambda
integer :: i,j
do j=1,N_states
do i=1,N_det
psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j)
enddo
call normalize(psi_coef(1,j), N_det)
enddo
SOFT_TOUCH psi_coef
end
logical function is_generable(det1, det2, Nint)

101
plugins/MRCC_Utils/multi_state.irp.f
View File

@ -1,101 +0,0 @@
subroutine multi_state(CI_electronic_energy_dressed_,CI_eigenvectors_dressed_,LDA)
implicit none
BEGIN_DOC
! Multi-state mixing
END_DOC
integer, intent(in) :: LDA
double precision, intent(inout) :: CI_electronic_energy_dressed_(N_states)
double precision, intent(inout) :: CI_eigenvectors_dressed_(LDA,N_states)
double precision, allocatable :: h(:,:,:), s(:,:), Psi(:,:), H_Psi(:,:,:), H_jj(:)
allocate( h(N_states,N_states,0:N_states), s(N_states,N_states) )
allocate( Psi(LDA,N_states), H_Psi(LDA,N_states,0:N_states) )
allocate (H_jj(LDA) )
! e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
integer :: i,j,k,istate
double precision :: U(N_states,N_states), Vt(N_states,N_states), D(N_states)
double precision, external :: diag_H_mat_elem
do istate=1,N_states
do i=1,N_det
H_jj(i) = diag_H_mat_elem(psi_det(1,1,i),N_int)
enddo
do i=1,N_det_ref
H_jj(idx_ref(i)) += delta_ii(istate,i)
enddo
do k=1,N_states
do i=1,N_det
Psi(i,k) = CI_eigenvectors_dressed_(i,k)
enddo
enddo
call H_u_0_mrcc_nstates(H_Psi(1,1,istate),Psi,H_jj,N_det,psi_det,N_int,istate,N_states,LDA)
do k=1,N_states
do i=1,N_states
double precision, external :: u_dot_v
h(i,k,istate) = u_dot_v(Psi(1,i), H_Psi(1,k,istate), N_det)
enddo
enddo
enddo
do k=1,N_states
do i=1,N_states
s(i,k) = u_dot_v(Psi(1,i), Psi(1,k), N_det)
enddo
enddo
print *, s(:,:)
print *, ''
h(:,:,0) = h(:,:,1)
do istate=2,N_states
U(:,:) = h(:,:,0)
call dgemm('N','N',N_states,N_states,N_states,1.d0,&
U, size(U,1), h(1,1,istate), size(h,1), 0.d0, &
h(1,1,0), size(Vt,1))
enddo
call svd(h(1,1,0), size(h,1), U, size(U,1), D, Vt, size(Vt,1), N_states, N_states)
do k=1,N_states
D(k) = D(k)**(1./dble(N_states))
if (D(k) > 0.d0) then
D(k) = -D(k)
endif
enddo
do j=1,N_states
do i=1,N_states
h(i,j,0) = 0.d0
do k=1,N_states
h(i,j,0) += U(i,k) * D(k) * Vt(k,j)
enddo
enddo
enddo
print *, h(:,:,0)
print *,''
integer :: LWORK, INFO
double precision, allocatable :: WORK(:)
LWORK=3*N_states
allocate (WORK(LWORK))
call dsygv(1, 'V', 'U', N_states, h(1,1,0), size(h,1), s, size(s,1), D, WORK, LWORK, INFO)
deallocate(WORK)
do j=1,N_states
do i=1,N_det
CI_eigenvectors_dressed_(i,j) = 0.d0
do k=1,N_states
CI_eigenvectors_dressed_(i,j) += Psi(i,k) * h(k,j,0)
enddo
enddo
CI_electronic_energy_dressed_(j) = D(j)
enddo
deallocate (h,s, H_jj)
deallocate( Psi, H_Psi )
end

2
plugins/Properties/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants Davidson
Determinants DavidsonUndressed

2
plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Psiref_Utils Davidson
Psiref_Utils DavidsonUndressed

1
plugins/UndressedMethod/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@

8
plugins/GPI2/README.rst → plugins/UndressedMethod/README.rst
View File

@ -1,8 +1,8 @@
=====
GASPI
=====
===============
UndressedMethod
===============
Providers for GASPI programs (with the GPI2 library).
Defines a null dressing vector
Needed Modules
==============

10
plugins/UndressedMethod/null_dressing_vector.irp.f Normal file
View File

@ -0,0 +1,10 @@
BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
&BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
implicit none
BEGIN_DOC
! Null dressing vectors
END_DOC
dressing_column_h(:,:) = 0.d0
dressing_column_s(:,:) = 0.d0
END_PROVIDER

2
plugins/mrcepa0/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ
Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ

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

@ -74,10 +74,8 @@ BEGIN_PROVIDER [ double precision, mrcc_norm_acc, (0:N_det_non_ref, N_states) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_sto, (N_states, N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_sto, (N_states, N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
integer :: gen, h, p, n, t, i, j, h1, h2, p1, p2, s1, s2, iproc
@ -94,10 +92,8 @@ END_PROVIDER
read(*,*) n_in_teeth
!n_in_teeth = 2
in_teeth_step = 1d0 / dfloat(n_in_teeth)
!double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref,N_det_ref) ]
!double precision :: delta_ii_mrcc_tmp, (N_states,N_det_ref) ]
!double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref,N_det_ref)
!double precision :: delta_ii_s2_mrcc_tmp(N_states, N_det_ref)
!double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref)
!double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref)
coefs = 0d0
coefs(:mrcc_teeth(1,1)-1) = 1d0
@ -144,15 +140,13 @@ END_PROVIDER
delta_ij_mrcc_sto = 0d0
delta_ii_mrcc_sto = 0d0
delta_ij_s2_mrcc_sto = 0d0
delta_ii_s2_mrcc_sto = 0d0
PROVIDE dij
provide hh_shortcut psi_det_size! lambda_mrcc
!$OMP PARALLEL DO default(none) schedule(dynamic) &
!$OMP shared(psi_ref, psi_non_ref, hh_exists, pp_exists, N_int, hh_shortcut) &
!$OMP shared(N_det_generators, coefs,N_det_non_ref, N_det_ref, delta_ii_mrcc_sto, delta_ij_mrcc_sto) &
!$OMP shared(contrib,psi_det_generators, delta_ii_s2_mrcc_sto, delta_ij_s2_mrcc_sto) &
!$OMP shared(N_det_generators, coefs,N_det_non_ref, delta_ij_mrcc_sto) &
!$OMP shared(contrib,psi_det_generators, delta_ij_s2_mrcc_sto) &
!$OMP private(i,j,curnorm,myCoef, h, n, mask, omask, buf, ok, iproc)
do gen= 1,N_det_generators
if(coefs(gen) == 0d0) cycle
@ -174,8 +168,8 @@ END_PROVIDER
end do
n = n - 1
if(n /= 0) then
call mrcc_part_dress(delta_ij_mrcc_sto, delta_ii_mrcc_sto, delta_ij_s2_mrcc_sto, &
delta_ii_s2_mrcc_sto, gen,n,buf,N_int,omask,myCoef,contrib)
call mrcc_part_dress(delta_ij_mrcc_sto, delta_ij_s2_mrcc_sto, &
gen,n,buf,N_int,omask,myCoef,contrib)
endif
end do
deallocate(buf)
@ -185,21 +179,17 @@ END_PROVIDER
curnorm = 0d0
do i=1,N_det_ref
do j=1,N_det_non_ref
curnorm += delta_ij_mrcc_sto(1, j, i)**2
curnorm += delta_ij_mrcc_sto(1,j)*delta_ij_mrcc_sto(1,j)
end do
end do
print *, "NORM DELTA ", curnorm**0.5d0
print *, "NORM DELTA ", dsqrt(curnorm)
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_cancel, (N_states, N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_cancel, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2_cancel, (N_states, N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_cancel, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
@ -216,15 +206,19 @@ END_PROVIDER
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2),inac, virt
integer, external :: get_index_in_psi_det_sorted_bit, searchDet,detCmp
logical, external :: is_in_wavefunction
double precision :: c0(N_states)
provide dij
delta_ij_cancel = 0d0
delta_ii_cancel = 0d0
do i_state = 1, N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
do i=1,N_det_ref
!$OMP PARALLEL DO default(shared) private(kk, k, blok, exc_Ik,det_tmp2,ok,deg,phase_Ik, l,ll) &
!$OMP private(contrib, contrib_s2, i_state)
!$OMP private(contrib, contrib_s2, i_state, c0)
do kk = 1, nlink(i)
k = det_cepa0_idx(linked(kk, i))
blok = blokMwen(kk, i)
@ -244,21 +238,10 @@ END_PROVIDER
do i_state = 1, N_states
contrib = (dij(j, l, i_state) - dij(i, k, i_state)) * delta_cas(i,j,i_state)! * Hla *phase_ia * phase_ik
contrib_s2 = dij(j, l, i_state) - dij(i, k, i_state)! * Sla*phase_ia * phase_ik
if(dabs(psi_ref_coef(i,i_state)).ge.1.d-3) then
!$OMP ATOMIC
delta_ij_cancel(i_state,l,i) += contrib
!$OMP ATOMIC
delta_ij_s2_cancel(i_state,l,i) += contrib_s2
!$OMP ATOMIC
delta_ii_cancel(i_state,i) -= contrib / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state)
!$OMP ATOMIC
delta_ii_s2_cancel(i_state,i) -= contrib_s2 / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state)
else
!$OMP ATOMIC
delta_ij_cancel(i_state,l,i) += contrib * 0.5d0
!$OMP ATOMIC
delta_ij_s2_cancel(i_state,l,i) += contrib_s2 * 0.5d0
endif
!$OMP ATOMIC
delta_ij_cancel(i_state,l) += contrib * psi_ref_coef(i,i_state) * c0(i_state)
!$OMP ATOMIC
delta_ij_s2_cancel(i_state,l) += contrib_s2* psi_ref_coef(i,i_state) * c0(i_state)
end do
end do
end do
@ -268,10 +251,8 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc, (N_states, N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc
@ -286,14 +267,12 @@ END_PROVIDER
contrib = 0d0
delta_ij_mrcc = 0d0
delta_ii_mrcc = 0d0
delta_ij_s2_mrcc = 0d0
delta_ii_s2_mrcc = 0d0
!$OMP PARALLEL DO default(none) schedule(dynamic) &
!$OMP shared(contrib,psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) &
!$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc, delta_ii_s2_mrcc, delta_ij_s2_mrcc) &
!$OMP shared(N_det_non_ref, N_det_ref, delta_ij_mrcc, delta_ij_s2_mrcc) &
!$OMP private(h, n, mask, omask, buf, ok, iproc)
do gen= 1, N_det_generators
allocate(buf(N_int, 2, N_det_non_ref))
@ -313,7 +292,7 @@ END_PROVIDER
n = n - 1
if(n /= 0) then
call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc, delta_ij_s2_mrcc, delta_ii_s2_mrcc, gen,n,buf,N_int,omask,1d0,contrib)
call mrcc_part_dress(delta_ij_mrcc, delta_ij_s2_mrcc, gen,n,buf,N_int,omask,1d0,contrib)
endif
end do
@ -324,20 +303,18 @@ END_PROVIDER
! subroutine blit(b1, b2)
! double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref)
! double precision :: b1(N_states,N_det_non_ref), b2(N_states,N_det_non_ref)
! b1 = b1 + b2
! end subroutine
subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,coef,contrib)
subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_buffer,Nint,key_mask,coef,contrib)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint
double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref)
double precision, intent(inout) :: delta_ii_(N_states,N_det_ref)
double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref)
double precision, intent(inout) :: delta_ii_s2_(N_states,N_det_ref)
double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref)
double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l,m
@ -399,6 +376,11 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
deallocate(microlist, idx_microlist)
double precision :: c0(N_states)
do i_state=1,N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
allocate (dIa_hla(N_states,N_det_non_ref), dIa_sla(N_states,N_det_non_ref))
! |I>
@ -436,8 +418,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
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(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)
@ -445,7 +427,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
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)
@ -467,9 +449,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
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))
call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_cache(k_sd))
!if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd)
!if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd)
enddo
! |I>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
@ -484,12 +466,12 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
! <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
! <I| /k\ |alpha>
! |l> = Exc(k -> alpha) |I>
@ -499,7 +481,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
logical :: ok
logical :: ok
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
do i_state=1,N_states
@ -510,7 +492,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
if (ok) then
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)
@ -522,40 +504,40 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
exit
endif
enddo
else if (perturbative_triples) then
! Linked
hka = hij_cache(idx_alpha(k_sd))
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
! Linked
hka = hij_cache(idx_alpha(k_sd))
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif
if (perturbative_triples.and. (degree2 == 1) ) then
call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
hka = hij_cache(idx_alpha(k_sd)) - hka
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
hka = hij_cache(idx_alpha(k_sd)) - hka
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif
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) = psi_ref_coef_inv(i_I,i_state)
enddo
@ -569,39 +551,17 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
enddo
enddo
do i_state=1,N_states
if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
p1 = 1
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
!$OMP ATOMIC
contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd,p1) += hdress
!$OMP ATOMIC
!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)
delta_ii_(i_state,p1) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd,p1) += sdress
!$OMP ATOMIC
!delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
delta_ii_s2_(i_state,p1) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
enddo
else
!stop "dress with coef < 1d-3"
delta_ii_(i_state,1) = 0.d0
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
p1 = 1
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd,p1) = delta_ij_(i_state,k_sd,p1) + 0.5d0*hdress
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd,p1) = delta_ij_s2_(i_state,k_sd,p1) + 0.5d0*sdress
enddo
endif
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
!$OMP ATOMIC
contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd) += hdress
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd) += sdress
enddo
enddo
enddo
enddo
@ -611,15 +571,13 @@ end
subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,contrib)
subroutine mrcc_part_dress_1c(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_buffer,Nint,key_mask,contrib)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint
double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref)
double precision, intent(inout) :: delta_ii_(N_states)
double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref)
double precision, intent(inout) :: delta_ii_s2_(N_states)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l,m
@ -715,6 +673,11 @@ subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_
end if
end if
double precision :: c0(N_states)
do i_state=1,N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
do i_alpha=1,N_tq
if(key_mask(1,1) /= 0) then
@ -850,39 +813,17 @@ subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_
enddo
enddo
do i_state=1,N_states
if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
p1 = 1
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
!$OMP ATOMIC
contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state)
contrib(i_state) += hdress * psi_ref_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd) += hdress
!$OMP ATOMIC
!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)
delta_ii_(i_state) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd) += sdress
!$OMP ATOMIC
!delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
delta_ii_s2_(i_state) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
enddo
else
!stop "dress with coef < 1d-3"
delta_ii_(i_state) = 0.d0
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
p1 = 1
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd) = delta_ij_(i_state,k_sd) + 0.5d0*hdress
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd) = delta_ij_s2_(i_state,k_sd) + 0.5d0*sdress
enddo
endif
enddo
enddo
enddo
@ -900,10 +841,8 @@ end
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_zmq, (N_states, N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_zmq, (N_states, N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij_mrcc_zmq, (N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_zmq, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
@ -917,9 +856,7 @@ end
delta_ij_mrcc_zmq = 0d0
delta_ii_mrcc_zmq = 0d0
delta_ij_s2_mrcc_zmq = 0d0
delta_ii_s2_mrcc_zmq = 0d0
!call random_seed()
E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion
@ -933,142 +870,67 @@ end
call ZMQ_mrcc(E_CI_before, mrcc, delta_ij_mrcc_zmq, delta_ij_s2_mrcc_zmq, abs(target_error))
mrcc_previous_E(:) = mrcc_E0_denominator(:)
do i=N_det_non_ref,1,-1
delta_ii_mrcc_zmq(:,1) -= delta_ij_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1)
delta_ii_s2_mrcc_zmq(:,1) -= delta_ij_s2_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1)
end do
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2, (N_states, N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
integer :: i, j, i_state
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc, 4=stoch
if(mrmode == 4) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_mrcc_sto(i_state,i)
delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_sto(i_state,i)
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_mrcc_sto(i_state,j,i)
delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_sto(i_state,j,i)
delta_ij(i_state,j) = delta_ij_mrcc_sto(i_state,j)
delta_ij_s2(i_state,j) = delta_ij_s2_mrcc_sto(i_state,j)
enddo
end do
end do
! else if(mrmode == 10) then
! do i = 1, N_det_ref
! do i_state = 1, N_states
! delta_ii(i_state,i)= delta_ii_mrsc2(i_state,i)
! delta_ii_s2(i_state,i)= delta_ii_s2_mrsc2(i_state,i)
! enddo
! do j = 1, N_det_non_ref
! do i_state = 1, N_states
! delta_ij(i_state,j,i) = delta_ij_mrsc2(i_state,j,i)
! delta_ij_s2(i_state,j,i) = delta_ij_s2_mrsc2(i_state,j,i)
! delta_ij(i_state,j) = delta_ij_mrsc2(i_state,j)
! delta_ij_s2(i_state,j) = delta_ij_s2_mrsc2(i_state,j)
! enddo
! end do
! end do
else if(mrmode == 5) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_mrcc_zmq(i_state,i)
delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_zmq(i_state,i)
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_mrcc_zmq(i_state,j,i)
delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_zmq(i_state,j,i)
delta_ij(i_state,j) = delta_ij_mrcc_zmq(i_state,j)
delta_ij_s2(i_state,j) = delta_ij_s2_mrcc_zmq(i_state,j)
enddo
end do
end do
else if(mrmode == 3) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
delta_ii_s2(i_state,i)= delta_ii_s2_mrcc(i_state,i)
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i)
delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc(i_state,j,i)
delta_ij(i_state,j) = delta_ij_mrcc(i_state,j)
delta_ij_s2(i_state,j) = delta_ij_s2_mrcc(i_state,j)
enddo
end do
end do
! =-=-= BEGIN STATE AVERAGE
! do i = 1, N_det_ref
! delta_ii(:,i)= delta_ii_mrcc(1,i)
! delta_ii_s2(:,i)= delta_ii_s2_mrcc(1,i)
! do i_state = 2, N_states
! delta_ii(:,i) += delta_ii_mrcc(i_state,i)
! delta_ii_s2(:,i) += delta_ii_s2_mrcc(i_state,i)
! enddo
! do j = 1, N_det_non_ref
! delta_ij(:,j,i) = delta_ij_mrcc(1,j,i)
! delta_ij_s2(:,j,i) = delta_ij_s2_mrcc(1,j,i)
! do i_state = 2, N_states
! delta_ij(:,j,i) += delta_ij_mrcc(i_state,j,i)
! delta_ij_s2(:,j,i) += delta_ij_s2_mrcc(i_state,j,i)
! enddo
! end do
! end do
! delta_ij = delta_ij * (1.d0/dble(N_states))
! delta_ii = delta_ii * (1.d0/dble(N_states))
! =-=-= END STATE AVERAGE
!
! do i = 1, N_det_ref
! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
! do j = 1, N_det_non_ref
! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
! end do
! end do
else if(mrmode == 2) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_old(i_state,i)
delta_ii_s2(i_state,i)= delta_ii_s2_old(i_state,i)
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
delta_ij_s2(i_state,j,i) = delta_ij_s2_old(i_state,j,i)
delta_ij(i_state,j) = delta_ij_old(i_state,j)
delta_ij_s2(i_state,j) = delta_ij_s2_old(i_state,j)
enddo
end do
end do
else if(mrmode == 1) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state)
delta_ii_s2(i_state,i)= delta_mrcepa0_ii_s2(i,i_state)
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state)
delta_ij_s2(i_state,j,i) = delta_mrcepa0_ij_s2(i,j,i_state)
delta_ij(i_state,j) = delta_mrcepa0_ij(j,i_state)
delta_ij_s2(i_state,j) = delta_mrcepa0_ij_s2(j,i_state)
enddo
end do
end do
else
stop "invalid mrmode"
end if
!if(mrmode == 2 .or. mrmode == 3) then
! do i = 1, N_det_ref
! do i_state = 1, N_states
! delta_ii(i_state,i) += delta_ii_cancel(i_state,i)
! enddo
! do j = 1, N_det_non_ref
! do i_state = 1, N_states
! delta_ij(i_state,j,i) += delta_ij_cancel(i_state,j,i)
! delta_ij(i_state,j) += delta_ij_cancel(i_state,j)
! enddo
! end do
! end do
!end if
END_PROVIDER
@ -1350,10 +1212,8 @@ subroutine getHP(a,h,p,Nint)
end subroutine
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_s2, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_s2, (N_det_ref,N_states) ]
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_s2, (N_det_non_ref,N_states) ]
use bitmasks
implicit none
@ -1361,7 +1221,7 @@ end subroutine
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref)
logical :: ok
double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
double precision :: contrib, contrib2, contrib_s2, contrib2_s2, HIIi, HJk, wall
double precision :: contrib, contrib_s2, HIIi, HJk, wall
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2)
integer(bit_kind),allocatable :: sortRef(:,:,:)
@ -1383,20 +1243,23 @@ end subroutine
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
enddo
double precision :: c0(N_states)
do i_state=1,N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
! To provide everything
contrib = dij(1, 1, 1)
delta_mrcepa0_ii(:,:) = 0d0
delta_mrcepa0_ij(:,:,:) = 0d0
delta_mrcepa0_ii_s2(:,:) = 0d0
delta_mrcepa0_ij_s2(:,:,:) = 0d0
delta_mrcepa0_ij(:,:) = 0d0
delta_mrcepa0_ij_s2(:,:) = 0d0
do i_state = 1, N_states
!$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) &
!$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) &
!$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ij_s2) &
!$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib_s2) &
!$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) &
!$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas, delta_cas_s2) &
!$OMP shared(notf,i_state, sortRef, sortRefIdx, dij)
!$OMP shared(notf,i_state, sortRef, sortRefIdx, dij,c0)
do blok=1,cepa0_shortcut(0)
do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
do II=1,N_det_ref
@ -1436,23 +1299,12 @@ end subroutine
!$OMP ATOMIC
notf = notf+1
! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
contrib = delta_cas(II, J, i_state)* dij(J, det_cepa0_idx(k), i_state)
contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then
contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
contrib2_s2 = contrib_s2 / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
!$OMP ATOMIC
delta_mrcepa0_ii(J,i_state) -= contrib2
delta_mrcepa0_ii_s2(J,i_state) -= contrib2_s2
else
contrib = contrib * 0.5d0
contrib_s2 = contrib_s2 * 0.5d0
end if
!$OMP ATOMIC
delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
delta_mrcepa0_ij_s2(J, det_cepa0_idx(i), i_state) += contrib_s2
delta_mrcepa0_ij(det_cepa0_idx(i), i_state) += contrib * c0(i_state) * psi_ref_coef(J,i_state)
delta_mrcepa0_ij_s2(det_cepa0_idx(i), i_state) += contrib_s2 * c0(i_state) * psi_ref_coef(J,i_state)
end do kloop
end do
@ -1467,8 +1319,7 @@ end subroutine
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ]
BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_non_ref,N_states) ]
use bitmasks
implicit none
@ -1476,7 +1327,7 @@ END_PROVIDER
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
logical :: ok
double precision :: phase_Ji, phase_Ik, phase_Ii
double precision :: contrib, contrib2, delta_IJk, HJk, HIk, HIl
double precision :: contrib, delta_IJk, HJk, HIk, HIl
integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
integer, allocatable :: idx_sorted_bit(:)
@ -1490,21 +1341,27 @@ END_PROVIDER
do i=1,N_det_non_ref
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
enddo
double precision :: c0(N_states)
do i_state=1,N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
do i_state = 1, N_states
delta_sub_ij(:,:,:) = 0d0
delta_sub_ii(:,:) = 0d0
delta_sub_ij(:,:) = 0d0
provide mo_bielec_integrals_in_map
!$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) &
!$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij) &
!$OMP private(i, J, k, degree, degree2, l, deg, ni) &
!$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) &
!$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
!$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
!$OMP private(det_tmp, det_tmp2, II, blok) &
!$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) &
!$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb)
!$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb,c0)
do i=1,N_det_non_ref
if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref
do J=1,N_det_ref
@ -1551,15 +1408,8 @@ END_PROVIDER
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
if(ok) cycle
contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)
if(dabs(psi_ref_coef(II,i_state)).ge.1.d-3) then
contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
!$OMP ATOMIC
delta_sub_ii(II,i_state) -= contrib2
else
contrib = contrib * 0.5d0
endif
!$OMP ATOMIC
delta_sub_ij(II, i, i_state) += contrib
delta_sub_ij(i, i_state) += contrib* c0(i_state) * psi_ref_coef(II,i_state)
end do
end do
end do

52
plugins/mrcepa0/dressing_slave.irp.f
View File

@ -402,17 +402,15 @@ end
subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii_s2_,delta_ij_s2_)
subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ij_,delta_ij_s2_)
use f77_zmq
implicit none
BEGIN_DOC
! Collects results from the AO integral calculation
END_DOC
double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref)
double precision,intent(inout) :: delta_ii_(N_states,N_det_ref)
double precision,intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref)
double precision,intent(inout) :: delta_ii_s2_(N_states,N_det_ref)
double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref)
double precision,intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref)
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
! integer :: j,l
@ -431,15 +429,18 @@ subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii
integer :: I_i, J, l, i_state, n(2), kk
integer,allocatable :: idx(:,:)
delta_ii_(:,:) = 0d0
delta_ij_(:,:,:) = 0d0
delta_ii_s2_(:,:) = 0d0
delta_ij_s2_(:,:,:) = 0d0
delta_ij_(:,:) = 0d0
delta_ij_s2_(:,:) = 0d0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
allocate ( delta(N_states,0:N_det_non_ref,2), delta_s2(N_states,0:N_det_non_ref,2) )
double precision :: c0(N_states)
do i_state=1,N_states
c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
enddo
allocate(idx(N_det_non_ref,2))
more = 1
do while (more == 1)
@ -449,34 +450,19 @@ subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii
do l=1, n(1)
do i_state=1,N_states
delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1)
delta_ij_s2_(i_state,idx(l,1),i_I) += delta_s2(i_state,l,1)
delta_ij_(i_state,idx(l,1)) += delta(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state)
delta_ij_s2_(i_state,idx(l,1)) += delta_s2(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state)
end do
end do
do l=1, n(2)
do i_state=1,N_states
delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2)
delta_ij_s2_(i_state,idx(l,2),J) += delta_s2(i_state,l,2)
delta_ij_(i_state,idx(l,2)) += delta(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state)
delta_ij_s2_(i_state,idx(l,2)) += delta_s2(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state)
end do
end do
if(n(1) /= 0) then
do i_state=1,N_states
delta_ii_(i_state,i_I) += delta(i_state,0,1)
delta_ii_s2_(i_state,i_I) += delta_s2(i_state,0,1)
end do
end if
if(n(2) /= 0) then
do i_state=1,N_states
delta_ii_(i_state,J) += delta(i_state,0,2)
delta_ii_s2_(i_state,J) += delta_s2(i_state,0,2)
end do
end if
if (task_id /= 0) then
integer, external :: zmq_delete_task
if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) == -1) then
@ -495,10 +481,8 @@ end
BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_s2_old, (N_states,N_det_ref) ]
BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref) ]
implicit none
integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2
@ -612,11 +596,11 @@ end
print *, nzer, ntot, float(nzer) / float(ntot)
provide nproc
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old,zmq_socket_pull)&
!$OMP SHARED(delta_ij_old,delta_ij_s2_old,zmq_socket_pull)&
!$OMP PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call mrsc2_dressing_collector(zmq_socket_pull,delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old)
call mrsc2_dressing_collector(zmq_socket_pull,delta_ij_old,delta_ij_s2_old)
else
call mrsc2_dressing_slave_inproc(i)
endif

4
plugins/mrcepa0/mrcc_stoch_routines.irp.f
View File

@ -35,6 +35,10 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error)
double precision :: w!(N_states)
integer, external :: add_task_to_taskserver
state_average_weight(:) = 0.d0
state_average_weight(mrcc_stoch_istate) = 1.d0
TOUCH state_average_weight
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral mrcc_weight psi_selectors

8
plugins/mrcepa0/mrcepa0_general.irp.f
View File

@ -14,8 +14,6 @@ subroutine run(N_st,energy)
integer :: n_it_mrcc_max
double precision :: thresh_mrcc
double precision, allocatable :: lambda(:)
allocate (lambda(N_states))
thresh_mrcc = thresh_dressed_ci
n_it_mrcc_max = n_it_max_dressed_ci
@ -34,7 +32,6 @@ subroutine run(N_st,energy)
E_new = 0.d0
delta_E = 1.d0
iteration = 0
lambda = 1.d0
do while (delta_E > thresh_mrcc)
iteration += 1
print *, '==============================================='
@ -45,12 +42,9 @@ subroutine run(N_st,energy)
do i=1,N_st
call write_double(6,ci_energy_dressed(i),"Energy")
enddo
call diagonalize_ci_dressed(lambda)
call diagonalize_ci_dressed
E_new = mrcc_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
! if (.true.) then
! provide delta_ij_mrcc_pouet
! endif
delta_E = (E_new - E_old)/dble(N_states)
print *, ''
call write_double(6,thresh_mrcc,"thresh_mrcc")

13
plugins/mrcepa0/run_mrcc_slave.irp.f
View File

@ -35,19 +35,16 @@ subroutine run_mrcc_slave(thread,iproc,energy)
integer(bit_kind) :: mask(N_int,2), omask(N_int,2)
double precision,allocatable :: delta_ij_loc(:,:,:)
double precision,allocatable :: delta_ii_loc(:,:)
!double precision,allocatable :: delta_ij_s2_loc(:,:,:)
!double precision,allocatable :: delta_ii_s2_loc(:,:)
integer :: h,p,n
logical :: ok
double precision :: contrib(N_states)
allocate(delta_ij_loc(N_states,N_det_non_ref,2) &
,delta_ii_loc(N_states,2))! &
allocate(delta_ij_loc(N_states,N_det_non_ref,2) )
!,delta_ij_s2_loc(N_states,N_det_non_ref,N_det_ref) &
!,delta_ii_s2_loc(N_states, N_det_ref))
allocate(abuf(N_int, 2, N_det_non_ref))
allocate(abuf(N_int, 2, N_det_non_ref))
allocate(mrcc_detail(N_states, N_det_generators))
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
@ -81,9 +78,7 @@ subroutine run_mrcc_slave(thread,iproc,energy)
contrib = 0d0
i_generator = ind(i_i_generator)
delta_ij_loc = 0d0
delta_ii_loc = 0d0
!delta_ij_s2_loc = 0d0
!delta_ii_s2_loc = 0d0
!call select_connected(i_generator,energy,mrcc_detail(1, i_i_generator),buf,subset)
!!!!!!!!!!!!!!!!!!!!!!
@ -102,7 +97,7 @@ subroutine run_mrcc_slave(thread,iproc,energy)
n = n - 1
if(n /= 0) then
call mrcc_part_dress_1c(delta_ij_loc(1,1,1), delta_ii_loc(1,1), delta_ij_loc(1,1,2), delta_ii_loc(1,2), &
call mrcc_part_dress_1c(delta_ij_loc(1,1,1), delta_ij_loc(1,1,2), &
i_generator,n,abuf,N_int,omask,contrib)
endif
end do

2
plugins/read_integral/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Integrals_Monoelec Integrals_Bielec
Integrals_Monoelec Integrals_Bielec Hartree_Fock

14
plugins/read_integral/print_integrals_mo.irp.f
View File

@ -44,14 +44,12 @@ program print_integrals
do l=1,mo_tot_num
do k=1,mo_tot_num
do j=l,mo_tot_num
do i=k,mo_tot_num
!if (i>=j) then
double precision :: get_mo_bielec_integral
integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
if (dabs(integral) > mo_integrals_threshold) then
write (iunit,'(4(I6,X),F20.15)') i,j,k,l, integral
endif
!end if
do i=max(j,k),mo_tot_num
double precision :: get_mo_bielec_integral
integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
if (dabs(integral) > mo_integrals_threshold) then
write (iunit,'(4(I6,X),E25.15)') i,j,k,l, integral
endif
enddo
enddo
enddo

39
plugins/read_integral/read_integrals_mo.irp.f
View File

@ -5,8 +5,44 @@ program read_integrals
! - nuclear_mo
! - bielec_mo
END_DOC
integer :: iunit
integer :: getunitandopen
integer :: i,j,n
PROVIDE ezfio_filename
call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("None")
logical :: has
call ezfio_has_mo_basis_mo_tot_num(has)
if (.not.has) then
iunit = getunitandopen('nuclear_mo','r')
n=0
do
read (iunit,*,end=12) i
n = max(n,i)
enddo
12 continue
close(iunit)
call ezfio_set_mo_basis_mo_tot_num(n)
call ezfio_has_ao_basis_ao_num(has)
mo_label = "None"
if (has) then
call huckel_guess
else
call ezfio_set_ao_basis_ao_num(n)
double precision, allocatable :: X(:,:)
allocate (X(n,n))
X = 0.d0
do i=1,n
X(i,i) = 1.d0
enddo
call ezfio_set_mo_basis_mo_coef(X)
call save_mos
endif
endif
call run
end
@ -69,9 +105,10 @@ subroutine run
13 continue
close(iunit)
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_values,0.d0)
call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals)
call map_sort(mo_integrals_map)
call map_unique(mo_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")

2
scripts/ezfio_interface/qp_edit_template
View File

@ -256,7 +256,7 @@ let spec =
let command =
Command.basic
Command.basic_spec
~summary: "Quantum Package command"
~readme:(fun () ->
"

27
scripts/ezfio_interface/upgrade_1.0_2.0.sh
View File

@ -1,27 +0,0 @@
#!/bin/bash
# Convert a old ezfio file (with option.irp.f ezfio_default)
# into a new EZFIO.cfg type
# Hartree Fock
# Changin the case, don't know if is needed or not
mv $1/Hartree_Fock $1/hartree_fock 2> /dev/null
mv $1/hartree_Fock/thresh_SCF $1/hartree_fock/thresh_scf 2> /dev/null
# BiInts
mv $1/bi_integrals $1/bielect_integrals 2> /dev/null
if [ -f $1/bielect_integrals/read_ao_integrals ]; then
if [ `cat $1/bielect_integrals/read_ao_integrals` -eq "True" ]
then
echo "Read" > $1/bielect_integrals/disk_access_ao_integrals
elif [ `cat bielect_integrals/write_ao_integrals` -eq "True" ]
then
echo "Write" > $1/bielect_integrals/disk_access_ao_integrals
else
echo "None" > $1/bielect_integrals/disk_access_ao_integrals
fi
fi

19
scripts/qp_upgrade_ocaml.sh Executable file
View File

@ -0,0 +1,19 @@
#!/bin/bash
OCAML_VERSION="4.06.0"
PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 ZMQ ppx_sexp_conv ppx_deriving"
if [[ -z ${QP_ROOT} ]]
then
print "The QP_ROOT environment variable is not set."
print "Please reload the quantum_package.rc file."
exit -1
fi
cd $QP_ROOT/ocaml
opam update
opam switch ${OCAML_VERSION}
eval `opam config env`
opam install -y ${PACKAGES} || echo "Upgrade failed. You can try running
configure ; $0"

2
src/Davidson/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants
Determinants DavidsonDressed

2
src/Davidson/diagonalize_CI.irp.f
View File

@ -65,7 +65,7 @@ END_PROVIDER
call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, &
size(CI_eigenvectors,1),CI_electronic_energy, &
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,6)
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,0)
else if (diag_algorithm == "Lapack") then

1
src/DavidsonDressed/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@

14
src/DavidsonDressed/README.rst Normal file
View File

@ -0,0 +1,14 @@
===============
DavidsonDressed
===============
Davidson with single-column dressing
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

103
src/Davidson/diagonalization_hs2.irp.f → src/DavidsonDressed/diagonalization_hs2_dressed.irp.f
View File

@ -1,4 +1,17 @@
subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit)
BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ]
implicit none
BEGIN_DOC
! Index of the dressed columns
END_DOC
integer :: i
double precision :: tmp
integer, external :: idamax
do i=1,N_states
dressed_column_idx(i) = idamax(size(psi_coef,1), psi_coef(1,i), 1)
enddo
END_PROVIDER
subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,dressing_state)
use bitmasks
implicit none
BEGIN_DOC
@ -15,41 +28,45 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
!
! N_st : Number of eigenstates
!
! iunit : Unit number for the I/O
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag)
double precision, allocatable :: H_jj(:)
integer, intent(in) :: dressing_state
double precision, allocatable :: H_jj(:), S2_jj(:)
double precision :: diag_H_mat_elem, diag_S_mat_elem
double precision, external :: diag_H_mat_elem, diag_S_mat_elem
integer :: i
ASSERT (N_st > 0)
ASSERT (sze > 0)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze) )
allocate(H_jj(sze),S2_jj(sze))
H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(sze,H_jj, dets_in,Nint) &
!$OMP PRIVATE(i)
!$OMP DO SCHEDULE(static)
do i=1,sze
do i=2,sze
H_jj(i) = diag_H_mat_elem(dets_in(1,1,i),Nint)
enddo
!$OMP END DO
!$OMP END PARALLEL
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
deallocate (H_jj)
if (dressing_state > 0) then
H_jj(dressed_column_idx(dressing_state)) += dressing_column_h(dressed_column_idx(dressing_state),dressing_state)
endif
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
deallocate (H_jj,S2_jj)
end
subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
use bitmasks
implicit none
BEGIN_DOC
@ -72,15 +89,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
!
! N_st_diag : Number of states in which H is diagonalized. Assumed > sze
!
! iunit : Unit for the I/O
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(in) :: H_jj(sze)
integer, intent(in) :: dressing_state
double precision, intent(inout) :: s2_out(N_st_diag)
integer, intent(in) :: iunit
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag)
@ -88,7 +103,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
integer :: i,j,k,l,m
logical :: converged
double precision :: u_dot_v, u_dot_u
double precision, external :: u_dot_v, u_dot_u
integer :: k_pairs, kl
@ -101,7 +116,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
character*(16384) :: write_buffer
double precision :: to_print(3,N_st)
double precision :: cpu, wall
integer :: shift, shift2, itermax
integer :: shift, shift2, itermax, istate
double precision :: r1, r2
logical :: state_ok(N_st_diag*davidson_sze_max)
include 'constants.include.F'
@ -117,35 +132,35 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse
call write_time(iunit)
call write_time(6)
call wall_time(wall)
call cpu_time(cpu)
write(iunit,'(A)') ''
write(iunit,'(A)') 'Davidson Diagonalization'
write(iunit,'(A)') '------------------------'
write(iunit,'(A)') ''
call write_int(iunit,N_st,'Number of states')
call write_int(iunit,N_st_diag,'Number of states in diagonalization')
call write_int(iunit,sze,'Number of determinants')
write(6,'(A)') ''
write(6,'(A)') 'Davidson Diagonalization'
write(6,'(A)') '------------------------'
write(6,'(A)') ''
call write_int(6,N_st,'Number of states')
call write_int(6,N_st_diag,'Number of states in diagonalization')
call write_int(6,sze,'Number of determinants')
r1 = 8.d0*(3.d0*dble(sze*N_st_diag*itermax+5.d0*(N_st_diag*itermax)**2 &
+ 4.d0*(N_st_diag*itermax)+nproc*(4.d0*N_det_alpha_unique+2.d0*N_st_diag*sze)))/(1024.d0**3)
call write_double(iunit, r1, 'Memory(Gb)')
write(iunit,'(A)') ''
call write_double(6, r1, 'Memory(Gb)')
write(6,'(A)') ''
write_buffer = '====='
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ =========== ==========='
enddo
write(iunit,'(A)') write_buffer(1:6+41*N_states)
write(6,'(A)') write_buffer(1:6+41*N_states)
write_buffer = 'Iter'
do i=1,N_st
write_buffer = trim(write_buffer)//' Energy S^2 Residual '
enddo
write(iunit,'(A)') write_buffer(1:6+41*N_states)
write(6,'(A)') write_buffer(1:6+41*N_states)
write_buffer = '====='
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ =========== ==========='
enddo
write(iunit,'(A)') write_buffer(1:6+41*N_states)
write(6,'(A)') write_buffer(1:6+41*N_states)
allocate( &
@ -225,7 +240,21 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
call H_S2_u_0_nstates_openmp(W(1,shift+1),S(1,shift+1),U(1,shift+1),N_st_diag,sze)
endif
if (dressing_state > 0) then
do istate=1,N_st_diag
l = dressed_column_idx(dressing_state)
do i=1,sze
W(i,shift+istate) += dressing_column_h(i,dressing_state) * U(l,shift+istate)
S(i,shift+istate) += dressing_column_s(i,dressing_state) * U(l,shift+istate)
W(l,shift+istate) += dressing_column_h(i,dressing_state) * U(i,shift+istate)
S(l,shift+istate) += dressing_column_s(i,dressing_state) * U(i,shift+istate)
enddo
W(l,shift+istate) -= dressing_column_h(l,dressing_state) * U(l,shift+istate)
S(l,shift+istate) -= dressing_column_s(l,dressing_state) * U(l,shift+istate)
enddo
endif
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
@ -399,7 +428,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
endif
enddo
write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st)
write(6,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st)
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
do k=1,N_st
if (residual_norm(k) > 1.e8) then
@ -429,9 +458,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ =========== ==========='
enddo
write(iunit,'(A)') trim(write_buffer)
write(iunit,'(A)') ''
call write_time(iunit)
write(6,'(A)') trim(write_buffer)
write(6,'(A)') ''
call write_time(6)
deallocate ( &
W, residual_norm, &
@ -443,6 +472,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
)
end
subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
use bitmasks
implicit none

210
src/DavidsonDressed/diagonalize_CI.irp.f Normal file
View File

@ -0,0 +1,210 @@
BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(6)
do j=1,min(N_det,N_states_diag)
CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
enddo
do j=1,min(N_det,N_states)
write(st,'(I4)') j
call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st))
call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ]
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
implicit none
double precision :: ovrlp,u_dot_v
integer :: i_good_state
integer, allocatable :: index_good_state_array(:)
logical, allocatable :: good_state_array(:)
double precision, allocatable :: s2_values_tmp(:)
integer :: i_other_state
double precision, allocatable :: eigenvectors(:,:), eigenvectors_s2(:,:), eigenvalues(:)
integer :: i_state
double precision :: e_0
integer :: i,j,k,mrcc_state
double precision, allocatable :: s2_eigvalues(:)
double precision, allocatable :: e_array(:)
integer, allocatable :: iorder(:)
PROVIDE threshold_davidson nthreads_davidson
! Guess values for the "N_states" states of the CI_eigenvectors_dressed
do j=1,min(N_states,N_det)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = psi_coef(i,j)
enddo
enddo
do j=min(N_states,N_det)+1,N_states_diag
do i=1,N_det
CI_eigenvectors_dressed(i,j) = 0.d0
enddo
enddo
if (diag_algorithm == "Davidson") then
allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),&
eigenvectors_s2(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),&
eigenvalues(size(CI_electronic_energy_dressed,1)))
do j=1,min(N_states,N_det)
do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
do mrcc_state=1,N_states
do j=mrcc_state,min(N_states,N_det)
do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
call davidson_diag_HS2(psi_det,eigenvectors, eigenvectors_s2, &
size(eigenvectors,1), &
eigenvalues,N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,&
mrcc_state)
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo
do k=N_states+1,N_states_diag
CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
CI_electronic_energy_dressed(k) = eigenvalues(k)
enddo
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors_dressed,1))
deallocate (eigenvectors,eigenvalues)
else if (diag_algorithm == "Lapack") then
allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_dressed,size(H_matrix_dressed,1),N_det)
CI_electronic_energy_dressed(:) = 0.d0
if (s2_eig) then
i_state = 0
allocate (s2_eigvalues(N_det))
allocate(index_good_state_array(N_det),good_state_array(N_det))
good_state_array = .False.
call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int, &
N_det,size(eigenvectors,1))
do j=1,N_det
! Select at least n_states states with S^2 values closed to "expected_s2"
if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
i_state +=1
index_good_state_array(i_state) = j
good_state_array(j) = .True.
endif
if(i_state.eq.N_states) then
exit
endif
enddo
if(i_state .ne.0)then
! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j))
enddo
CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j))
CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j))
enddo
i_other_state = 0
do j = 1, N_det
if(good_state_array(j))cycle
i_other_state +=1
if(i_state+i_other_state.gt.n_states_diag)then
exit
endif
do i=1,N_det
CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j)
CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
enddo
else
print*,''
print*,'!!!!!!!! WARNING !!!!!!!!!'
print*,' Within the ',N_det,'determinants selected'
print*,' and the ',N_states_diag,'states requested'
print*,' We did not find any state with S^2 values close to ',expected_s2
print*,' We will then set the first N_states eigenvectors of the H matrix'
print*,' as the CI_eigenvectors_dressed'
print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
print*,''
do j=1,min(N_states_diag,N_det)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(j) = eigenvalues(j)
CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j)
enddo
endif
deallocate(index_good_state_array,good_state_array)
deallocate(s2_eigvalues)
else
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states_diag),size(eigenvectors,1))
! Select the "N_states_diag" states of lowest energy
do j=1,min(N_det,N_states_diag)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy_dressed(j) = eigenvalues(j)
enddo
endif
deallocate(eigenvectors,eigenvalues)
endif
END_PROVIDER
subroutine diagonalize_CI_dressed
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
enddo
enddo
SOFT_TOUCH psi_coef
end
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
implicit none
BEGIN_DOC
! Dressed H with Delta_ij
END_DOC
integer :: i, j,istate,ii,jj
do istate = 1,N_states
do j=1,N_det
do i=1,N_det
h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j)
enddo
enddo
i = dressed_column_idx(istate)
do j = 1, N_det
h_matrix_dressed(i,j,istate) += dressing_column_h(j,istate)
h_matrix_dressed(j,i,istate) += dressing_column_h(j,istate)
enddo
h_matrix_dressed(i,i,istate) -= dressing_column_h(i,istate)
enddo
END_PROVIDER

1
src/DavidsonUndressed/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Davidson UndressedMethod

14
src/DavidsonUndressed/README.rst Normal file
View File

@ -0,0 +1,14 @@
=================
DavidsonUndressed
=================
Module for main files with undressed Davidson
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

0
src/Davidson/davidson_slave.irp.f → src/DavidsonUndressed/davidson_slave.irp.f
View File

0
src/Davidson/diagonalize_restart_and_save_all_states.irp.f → src/DavidsonUndressed/diag_restart_save_all_states.irp.f
View File

0
src/Davidson/diagonalize_restart_and_save_lowest_state.irp.f → src/DavidsonUndressed/diag_restart_save_lowest_state.irp.f
View File

0
src/Davidson/diagonalize_restart_and_save_one_state.irp.f → src/DavidsonUndressed/diag_restart_save_one_state.irp.f
View File

0
src/Davidson/guess_lowest_state.irp.f → src/DavidsonUndressed/guess_lowest_state.irp.f
View File

0
src/Davidson/print_H_matrix_restart.irp.f → src/DavidsonUndressed/print_H_matrix_restart.irp.f
View File

0
src/Davidson/print_energy.irp.f → src/DavidsonUndressed/print_energy.irp.f
View File

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

@ -368,13 +368,13 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
END_DOC
logical :: exists
state_average_weight = 1.d0
state_average_weight(:) = 1.d0
call ezfio_has_determinants_state_average_weight(exists)
if (exists) then
call ezfio_get_determinants_state_average_weight(state_average_weight)
endif
state_average_weight = state_average_weight+1.d-31
state_average_weight = state_average_weight/(sum(state_average_weight))
state_average_weight(:) = state_average_weight(:)+1.d-31
state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:)))
END_PROVIDER

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

@ -225,34 +225,6 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
END_PROVIDER
subroutine update_psi_average_norm_contrib(w)
implicit none
BEGIN_DOC
! Compute psi_average_norm_contrib for different state average weights w(:)
END_DOC
double precision, intent(in) :: w(N_states)
double precision :: w0(N_states), f
w0(:) = w(:)/sum(w(:))
integer :: i,j,k
do i=1,N_det
psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*w(1)
enddo
do k=2,N_states
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + &
psi_coef(i,k)*psi_coef(i,k)*w(k)
enddo
enddo
f = 1.d0/sum(psi_average_norm_contrib(1:N_det))
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f
enddo
SOFT_TOUCH psi_average_norm_contrib
end subroutine
BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
implicit none
BEGIN_DOC
@ -260,14 +232,12 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
END_DOC
integer :: i,j,k
double precision :: f
f = 1.d0/dble(N_states)
do i=1,N_det
psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f
enddo
do k=2,N_states
psi_average_norm_contrib(:) = 0.d0
do k=1,N_states
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + &
psi_coef(i,k)*psi_coef(i,k)*f
psi_coef(i,k)*psi_coef(i,k)*state_average_weight(k)
enddo
enddo
f = 1.d0/sum(psi_average_norm_contrib(1:N_det))