mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-30 15:15:38 +01:00
Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev
This commit is contained in:
commit
1824e23ab4
@ -206,7 +206,7 @@ ZeroMQ and its Fortran binding
|
||||
|
||||
.. code:: bash
|
||||
|
||||
cp f77_zmq_free.h ${QP_ROOT}/src/ZMQ/f77_zmq.h
|
||||
cp f77_zmq_free.h ${QP_ROOT}/src/zmq/f77_zmq.h
|
||||
|
||||
|
||||
Zlib
|
||||
|
||||
63
config/ifort_2019_mpi_rome.cfg
Normal file
63
config/ifort_2019_mpi_rome.cfg
Normal file
@ -0,0 +1,63 @@
|
||||
# Common flags
|
||||
##############
|
||||
#
|
||||
# -mkl=[parallel|sequential] : Use the MKL library
|
||||
# --ninja : Allow the utilisation of ninja. It is mandatory !
|
||||
# --align=32 : Align all provided arrays on a 32-byte boundary
|
||||
#
|
||||
[COMMON]
|
||||
FC : mpiifort -fpic
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED
|
||||
|
||||
# Global options
|
||||
################
|
||||
#
|
||||
# 1 : Activate
|
||||
# 0 : Deactivate
|
||||
#
|
||||
[OPTION]
|
||||
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
CACHE : 0 ; Enable cache_compile.py
|
||||
OPENMP : 1 ; Append OpenMP flags
|
||||
|
||||
# Optimization flags
|
||||
####################
|
||||
#
|
||||
# -xHost : Compile a binary optimized for the current architecture
|
||||
# -O2 : O3 not better than O2.
|
||||
# -ip : Inter-procedural optimizations
|
||||
# -ftz : Flushes denormal results to zero
|
||||
#
|
||||
[OPT]
|
||||
FC : -traceback -shared-intel
|
||||
FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer
|
||||
|
||||
# Profiling flags
|
||||
#################
|
||||
#
|
||||
[PROFILE]
|
||||
FC : -p -g
|
||||
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
|
||||
|
||||
# Debugging flags
|
||||
#################
|
||||
#
|
||||
# -traceback : Activate backtrace on runtime
|
||||
# -fpe0 : All floating point exaceptions
|
||||
# -C : Checks uninitialized variables, array subscripts, etc...
|
||||
# -g : Extra debugging information
|
||||
# -xSSE2 : Valgrind needs a very simple x86 executable
|
||||
#
|
||||
[DEBUG]
|
||||
FC : -g -traceback
|
||||
FCFLAGS : -xSSE2 -C -fpe0 -implicitnone
|
||||
|
||||
# OpenMP flags
|
||||
#################
|
||||
#
|
||||
[OPENMP]
|
||||
FC : -qopenmp
|
||||
IRPF90_FLAGS : --openmp
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
#
|
||||
[COMMON]
|
||||
FC : ifort -fpic
|
||||
LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DINTEL
|
||||
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
#
|
||||
[COMMON]
|
||||
FC : mpiifort -fpic
|
||||
LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL
|
||||
|
||||
|
||||
63
config/ifort_2021_mpi_rome.cfg
Normal file
63
config/ifort_2021_mpi_rome.cfg
Normal file
@ -0,0 +1,63 @@
|
||||
# Common flags
|
||||
##############
|
||||
#
|
||||
# -mkl=[parallel|sequential] : Use the MKL library
|
||||
# --ninja : Allow the utilisation of ninja. It is mandatory !
|
||||
# --align=32 : Align all provided arrays on a 32-byte boundary
|
||||
#
|
||||
[COMMON]
|
||||
FC : mpiifort -fpic
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DINTEL
|
||||
|
||||
# Global options
|
||||
################
|
||||
#
|
||||
# 1 : Activate
|
||||
# 0 : Deactivate
|
||||
#
|
||||
[OPTION]
|
||||
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
CACHE : 0 ; Enable cache_compile.py
|
||||
OPENMP : 1 ; Append OpenMP flags
|
||||
|
||||
# Optimization flags
|
||||
####################
|
||||
#
|
||||
# -xHost : Compile a binary optimized for the current architecture
|
||||
# -O2 : O3 not better than O2.
|
||||
# -ip : Inter-procedural optimizations
|
||||
# -ftz : Flushes denormal results to zero
|
||||
#
|
||||
[OPT]
|
||||
FC : -traceback -shared-intel
|
||||
FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer
|
||||
|
||||
# Profiling flags
|
||||
#################
|
||||
#
|
||||
[PROFILE]
|
||||
FC : -p -g
|
||||
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
|
||||
|
||||
# Debugging flags
|
||||
#################
|
||||
#
|
||||
# -traceback : Activate backtrace on runtime
|
||||
# -fpe0 : All floating point exaceptions
|
||||
# -C : Checks uninitialized variables, array subscripts, etc...
|
||||
# -g : Extra debugging information
|
||||
# -xSSE2 : Valgrind needs a very simple x86 executable
|
||||
#
|
||||
[DEBUG]
|
||||
FC : -g -traceback
|
||||
FCFLAGS : -xSSE2 -C -fpe0 -implicitnone
|
||||
|
||||
# OpenMP flags
|
||||
#################
|
||||
#
|
||||
[OPENMP]
|
||||
FC : -qopenmp
|
||||
IRPF90_FLAGS : --openmp
|
||||
|
||||
63
config/ifort_2021_rome.cfg
Normal file
63
config/ifort_2021_rome.cfg
Normal file
@ -0,0 +1,63 @@
|
||||
# Common flags
|
||||
##############
|
||||
#
|
||||
# -mkl=[parallel|sequential] : Use the MKL library
|
||||
# --ninja : Allow the utilisation of ninja. It is mandatory !
|
||||
# --align=32 : Align all provided arrays on a 32-byte boundary
|
||||
#
|
||||
[COMMON]
|
||||
FC : ifort -fpic
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DINTEL
|
||||
|
||||
# Global options
|
||||
################
|
||||
#
|
||||
# 1 : Activate
|
||||
# 0 : Deactivate
|
||||
#
|
||||
[OPTION]
|
||||
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
CACHE : 0 ; Enable cache_compile.py
|
||||
OPENMP : 1 ; Append OpenMP flags
|
||||
|
||||
# Optimization flags
|
||||
####################
|
||||
#
|
||||
# -xHost : Compile a binary optimized for the current architecture
|
||||
# -O2 : O3 not better than O2.
|
||||
# -ip : Inter-procedural optimizations
|
||||
# -ftz : Flushes denormal results to zero
|
||||
#
|
||||
[OPT]
|
||||
FC : -traceback -shared-intel
|
||||
FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer
|
||||
|
||||
# Profiling flags
|
||||
#################
|
||||
#
|
||||
[PROFILE]
|
||||
FC : -p -g
|
||||
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
|
||||
|
||||
# Debugging flags
|
||||
#################
|
||||
#
|
||||
# -traceback : Activate backtrace on runtime
|
||||
# -fpe0 : All floating point exaceptions
|
||||
# -C : Checks uninitialized variables, array subscripts, etc...
|
||||
# -g : Extra debugging information
|
||||
# -xSSE2 : Valgrind needs a very simple x86 executable
|
||||
#
|
||||
[DEBUG]
|
||||
FC : -g -traceback
|
||||
FCFLAGS : -xSSE2 -C -fpe0 -implicitnone
|
||||
|
||||
# OpenMP flags
|
||||
#################
|
||||
#
|
||||
[OPENMP]
|
||||
FC : -qopenmp
|
||||
IRPF90_FLAGS : --openmp
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
#
|
||||
[COMMON]
|
||||
FC : ifort -fpic
|
||||
LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DINTEL
|
||||
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
#
|
||||
[COMMON]
|
||||
FC : mpiifort -fpic
|
||||
LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL
|
||||
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
#
|
||||
[COMMON]
|
||||
FC : ifort -fpic
|
||||
LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps
|
||||
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
|
||||
IRPF90 : irpf90
|
||||
IRPF90_FLAGS : --ninja --align=64 -DINTEL
|
||||
|
||||
|
||||
@ -277,6 +277,16 @@ let run ?o b au c d m p cart xyz_file =
|
||||
) nuclei
|
||||
in
|
||||
|
||||
let z_core =
|
||||
List.map (fun x ->
|
||||
Positive_int.to_int x.Pseudo.n_elec
|
||||
|> float_of_int
|
||||
) pseudo
|
||||
in
|
||||
let nucl_num = (List.length z_core) in
|
||||
Ezfio.set_pseudo_nucl_charge_remove (Ezfio.ezfio_array_of_list
|
||||
~rank:1 ~dim:[| nucl_num |] ~data:z_core);
|
||||
|
||||
let molecule =
|
||||
let n_elec_to_remove =
|
||||
List.fold_left (fun accu x ->
|
||||
|
||||
@ -12,21 +12,21 @@ double precision function ao_value(i,r)
|
||||
integer :: power_ao(3)
|
||||
double precision :: accu,dx,dy,dz,r2
|
||||
num_ao = ao_nucl(i)
|
||||
power_ao(1:3)= ao_power(i,1:3)
|
||||
center_ao(1:3) = nucl_coord(num_ao,1:3)
|
||||
dx = (r(1) - center_ao(1))
|
||||
dy = (r(2) - center_ao(2))
|
||||
dz = (r(3) - center_ao(3))
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
dx = dx**power_ao(1)
|
||||
dy = dy**power_ao(2)
|
||||
dz = dz**power_ao(3)
|
||||
! power_ao(1:3)= ao_power(i,1:3)
|
||||
! center_ao(1:3) = nucl_coord(num_ao,1:3)
|
||||
! dx = (r(1) - center_ao(1))
|
||||
! dy = (r(2) - center_ao(2))
|
||||
! dz = (r(3) - center_ao(3))
|
||||
! r2 = dx*dx + dy*dy + dz*dz
|
||||
! dx = dx**power_ao(1)
|
||||
! dy = dy**power_ao(2)
|
||||
! dz = dz**power_ao(3)
|
||||
|
||||
accu = 0.d0
|
||||
do m=1,ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
|
||||
enddo
|
||||
! do m=1,ao_prim_num(i)
|
||||
! beta = ao_expo_ordered_transp(m,i)
|
||||
! accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
|
||||
! enddo
|
||||
ao_value = accu * dx * dy * dz
|
||||
|
||||
end
|
||||
|
||||
@ -268,6 +268,21 @@ subroutine print_spindet(string,Nint)
|
||||
|
||||
end
|
||||
|
||||
subroutine print_det_one_dimension(string,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Subroutine to print the content of a determinant using the '+-' notation
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: string(Nint)
|
||||
character*(2048) :: output(1)
|
||||
|
||||
call bitstring_to_str( output(1), string, Nint )
|
||||
print *, trim(output(1))
|
||||
|
||||
end
|
||||
|
||||
logical function is_integer_in_string(bite,string,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
@ -1,9 +1,3 @@
|
||||
[pert_2rdm]
|
||||
type: logical
|
||||
doc: If true, computes the one- and two-body rdms with perturbation theory
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[save_wf_after_selection]
|
||||
type: logical
|
||||
doc: If true, saves the wave function after the selection, before the diagonalization
|
||||
|
||||
@ -2,5 +2,4 @@ perturbation
|
||||
zmq
|
||||
mpi
|
||||
iterations
|
||||
two_body_rdm
|
||||
csf
|
||||
|
||||
@ -70,8 +70,8 @@ subroutine run_cipsi
|
||||
|
||||
do while ( &
|
||||
(N_det < N_det_max) .and. &
|
||||
(maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. &
|
||||
(maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. &
|
||||
(sum(abs(pt2_data % pt2(1:N_states)) * state_average_weight(1:N_states)) > pt2_max) .and. &
|
||||
(sum(abs(pt2_data % variance(1:N_states)) * state_average_weight(1:N_states)) > variance_max) .and. &
|
||||
(correlation_energy_ratio <= correlation_energy_ratio_max) &
|
||||
)
|
||||
write(*,'(A)') '--------------------------------------------------------------------------------'
|
||||
|
||||
@ -1,183 +0,0 @@
|
||||
|
||||
use bitmasks
|
||||
use omp_lib
|
||||
|
||||
BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock]
|
||||
use f77_zmq
|
||||
implicit none
|
||||
call omp_init_lock(pert_2rdm_lock)
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, n_orb_pert_rdm]
|
||||
implicit none
|
||||
n_orb_pert_rdm = n_act_orb
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)]
|
||||
implicit none
|
||||
list_orb_reverse_pert_rdm = list_act_reverse
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
|
||||
implicit none
|
||||
list_orb_pert_rdm = list_act
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)]
|
||||
implicit none
|
||||
pert_2rdm_provider = 0.d0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection)
|
||||
use bitmasks
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n_det_connection
|
||||
double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection)
|
||||
integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection)
|
||||
integer, intent(in) :: i_generator, sp, h1, h2
|
||||
double precision, intent(in) :: mat(N_states, mo_num, mo_num)
|
||||
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num)
|
||||
double precision, intent(in) :: fock_diag_tmp(mo_num)
|
||||
double precision, intent(in) :: E0(N_states)
|
||||
type(pt2_type), intent(inout) :: pt2_data
|
||||
type(selection_buffer), intent(inout) :: buf
|
||||
logical :: ok
|
||||
integer :: s1, s2, p1, p2, ib, j, istate, jstate
|
||||
integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
|
||||
double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states)
|
||||
double precision, external :: diag_H_mat_elem_fock
|
||||
double precision :: E_shift
|
||||
|
||||
logical, external :: detEq
|
||||
double precision, allocatable :: values(:)
|
||||
integer, allocatable :: keys(:,:)
|
||||
integer :: nkeys
|
||||
integer :: sze_buff
|
||||
sze_buff = 5 * mo_num ** 2
|
||||
allocate(keys(4,sze_buff),values(sze_buff))
|
||||
nkeys = 0
|
||||
if(sp == 3) then
|
||||
s1 = 1
|
||||
s2 = 2
|
||||
else
|
||||
s1 = sp
|
||||
s2 = sp
|
||||
end if
|
||||
call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
|
||||
E_shift = 0.d0
|
||||
|
||||
if (h0_type == 'CFG') then
|
||||
j = det_to_configuration(i_generator)
|
||||
E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j)
|
||||
endif
|
||||
|
||||
do p1=1,mo_num
|
||||
if(bannedOrb(p1, s1)) cycle
|
||||
ib = 1
|
||||
if(sp /= 3) ib = p1+1
|
||||
|
||||
do p2=ib,mo_num
|
||||
|
||||
! -----
|
||||
! /!\ Generating only single excited determinants doesn't work because a
|
||||
! determinant generated by a single excitation may be doubly excited wrt
|
||||
! to a determinant of the future. In that case, the determinant will be
|
||||
! detected as already generated when generating in the future with a
|
||||
! double excitation.
|
||||
!
|
||||
! if (.not.do_singles) then
|
||||
! if ((h1 == p1) .or. (h2 == p2)) then
|
||||
! cycle
|
||||
! endif
|
||||
! endif
|
||||
!
|
||||
! if (.not.do_doubles) then
|
||||
! if ((h1 /= p1).and.(h2 /= p2)) then
|
||||
! cycle
|
||||
! endif
|
||||
! endif
|
||||
! -----
|
||||
|
||||
if(bannedOrb(p2, s2)) cycle
|
||||
if(banned(p1,p2)) cycle
|
||||
|
||||
|
||||
if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle
|
||||
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
|
||||
|
||||
if (do_only_cas) then
|
||||
integer, external :: number_of_holes, number_of_particles
|
||||
if (number_of_particles(det)>0) then
|
||||
cycle
|
||||
endif
|
||||
if (number_of_holes(det)>0) then
|
||||
cycle
|
||||
endif
|
||||
endif
|
||||
|
||||
if (do_ddci) then
|
||||
logical, external :: is_a_two_holes_two_particles
|
||||
if (is_a_two_holes_two_particles(det)) then
|
||||
cycle
|
||||
endif
|
||||
endif
|
||||
|
||||
if (do_only_1h1p) then
|
||||
logical, external :: is_a_1h1p
|
||||
if (.not.is_a_1h1p(det)) cycle
|
||||
endif
|
||||
|
||||
|
||||
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
|
||||
|
||||
sum_e_pert = 0d0
|
||||
integer :: degree
|
||||
call get_excitation_degree(det,HF_bitmask,degree,N_int)
|
||||
if(degree == 2)cycle
|
||||
do istate=1,N_states
|
||||
delta_E = E0(istate) - Hii + E_shift
|
||||
alpha_h_psi = mat(istate, p1, p2)
|
||||
val = alpha_h_psi + alpha_h_psi
|
||||
tmp = dsqrt(delta_E * delta_E + val * val)
|
||||
if (delta_E < 0.d0) then
|
||||
tmp = -tmp
|
||||
endif
|
||||
e_pert = 0.5d0 * (tmp - delta_E)
|
||||
coef(istate) = e_pert / alpha_h_psi
|
||||
print*,e_pert,coef,alpha_h_psi
|
||||
pt2_data % pt2(istate) += e_pert
|
||||
pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
|
||||
enddo
|
||||
|
||||
do istate=1,N_states
|
||||
alpha_h_psi = mat(istate, p1, p2)
|
||||
e_pert = coef(istate) * alpha_h_psi
|
||||
do jstate=1,N_states
|
||||
pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate)
|
||||
enddo
|
||||
|
||||
if (weight_selection /= 5) then
|
||||
! Energy selection
|
||||
sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
|
||||
|
||||
else
|
||||
! Variance selection
|
||||
sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
|
||||
endif
|
||||
end do
|
||||
call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
|
||||
|
||||
if(sum_e_pert <= buf%mini) then
|
||||
call add_to_selection_buffer(buf, det, sum_e_pert)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
|
||||
end
|
||||
|
||||
|
||||
@ -131,8 +131,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
|
||||
PROVIDE psi_det_hii selection_weight pseudo_sym
|
||||
PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
|
||||
PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max
|
||||
PROVIDE list_act list_inact list_core list_virt list_del seniority_max
|
||||
PROVIDE excitation_beta_max excitation_alpha_max excitation_max
|
||||
|
||||
if (h0_type == 'CFG') then
|
||||
PROVIDE psi_configuration_hii det_to_configuration
|
||||
|
||||
@ -464,14 +464,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
|
||||
allocate (fullminilist (N_int, 2, fullinteresting(0)), &
|
||||
minilist (N_int, 2, interesting(0)) )
|
||||
if(pert_2rdm)then
|
||||
allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
|
||||
do i=1,fullinteresting(0)
|
||||
do j = 1, N_states
|
||||
coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
! if(pert_2rdm)then
|
||||
! allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
|
||||
! do i=1,fullinteresting(0)
|
||||
! do j = 1, N_states
|
||||
! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
|
||||
! enddo
|
||||
! enddo
|
||||
! endif
|
||||
|
||||
do i=1,fullinteresting(0)
|
||||
do k=1,N_int
|
||||
@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
|
||||
call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
|
||||
|
||||
if(.not.pert_2rdm)then
|
||||
! if(.not.pert_2rdm)then
|
||||
call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
|
||||
else
|
||||
call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
|
||||
endif
|
||||
! else
|
||||
! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
|
||||
! endif
|
||||
end if
|
||||
enddo
|
||||
if(s1 /= s2) monoBdo = .false.
|
||||
enddo
|
||||
deallocate(fullminilist,minilist)
|
||||
if(pert_2rdm)then
|
||||
deallocate(coef_fullminilist_rev)
|
||||
endif
|
||||
! if(pert_2rdm)then
|
||||
! deallocate(coef_fullminilist_rev)
|
||||
! endif
|
||||
enddo
|
||||
enddo
|
||||
deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
|
||||
|
||||
@ -69,8 +69,8 @@ subroutine run_stochastic_cipsi
|
||||
|
||||
do while ( &
|
||||
(N_det < N_det_max) .and. &
|
||||
(maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. &
|
||||
(maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. &
|
||||
(sum(abs(pt2_data % pt2(1:N_states)) * state_average_weight(1:N_states)) > pt2_max) .and. &
|
||||
(sum(abs(pt2_data % variance(1:N_states)) * state_average_weight(1:N_states)) > variance_max) .and. &
|
||||
(correlation_energy_ratio <= correlation_energy_ratio_max) &
|
||||
)
|
||||
write(*,'(A)') '--------------------------------------------------------------------------------'
|
||||
|
||||
@ -1,223 +0,0 @@
|
||||
use bitmasks
|
||||
|
||||
subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
|
||||
implicit none
|
||||
integer, intent(in) :: n_det_connection,sze_buff
|
||||
double precision, intent(in) :: coef(N_states)
|
||||
integer(bit_kind), intent(in) :: det(N_int,2)
|
||||
integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection)
|
||||
double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection)
|
||||
integer, intent(inout) :: keys(4,sze_buff),nkeys
|
||||
double precision, intent(inout) :: values(sze_buff)
|
||||
integer :: i,j
|
||||
integer :: exc(0:2,2,2)
|
||||
integer :: degree
|
||||
double precision :: phase, contrib
|
||||
do i = 1, n_det_connection
|
||||
call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int)
|
||||
if(degree.gt.2)cycle
|
||||
contrib = 0.d0
|
||||
do j = 1, N_states
|
||||
contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j)
|
||||
enddo
|
||||
! case of single excitations
|
||||
if(degree == 1)then
|
||||
if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then
|
||||
call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
|
||||
nkeys = 0
|
||||
endif
|
||||
call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff)
|
||||
else
|
||||
!! case of double excitations
|
||||
! if (nkeys + 4 .ge. sze_buff)then
|
||||
! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
|
||||
! nkeys = 0
|
||||
! endif
|
||||
! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
|
||||
endif
|
||||
enddo
|
||||
!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
|
||||
!nkeys = 0
|
||||
|
||||
end
|
||||
|
||||
subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
|
||||
implicit none
|
||||
integer, intent(in) :: sze_buff
|
||||
integer(bit_kind), intent(in) :: det1(N_int,2)
|
||||
integer(bit_kind), intent(in) :: det2(N_int,2)
|
||||
integer,intent(in) :: exc(0:2,2,2)
|
||||
double precision,intent(in) :: phase, contrib
|
||||
integer, intent(inout) :: nkeys, keys(4,sze_buff)
|
||||
double precision, intent(inout):: values(sze_buff)
|
||||
|
||||
integer :: occ(N_int*bit_kind_size,2)
|
||||
integer :: n_occ_ab(2),ispin,other_spin
|
||||
integer :: h1,h2,p1,p2,i
|
||||
call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int)
|
||||
|
||||
if (exc(0,1,1) == 1) then
|
||||
! Mono alpha
|
||||
h1 = exc(1,1,1)
|
||||
p1 = exc(1,2,1)
|
||||
ispin = 1
|
||||
other_spin = 2
|
||||
else
|
||||
! Mono beta
|
||||
h1 = exc(1,1,2)
|
||||
p1 = exc(1,2,2)
|
||||
ispin = 2
|
||||
other_spin = 1
|
||||
endif
|
||||
if(list_orb_reverse_pert_rdm(h1).lt.0)return
|
||||
h1 = list_orb_reverse_pert_rdm(h1)
|
||||
if(list_orb_reverse_pert_rdm(p1).lt.0)return
|
||||
p1 = list_orb_reverse_pert_rdm(p1)
|
||||
!update the alpha/beta part
|
||||
do i = 1, n_occ_ab(other_spin)
|
||||
h2 = occ(i,other_spin)
|
||||
if(list_orb_reverse_pert_rdm(h2).lt.0)return
|
||||
h2 = list_orb_reverse_pert_rdm(h2)
|
||||
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h1
|
||||
keys(2,nkeys) = h2
|
||||
keys(3,nkeys) = p1
|
||||
keys(4,nkeys) = h2
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h2
|
||||
keys(2,nkeys) = h1
|
||||
keys(3,nkeys) = h2
|
||||
keys(4,nkeys) = p1
|
||||
enddo
|
||||
!update the same spin part
|
||||
!do i = 1, n_occ_ab(ispin)
|
||||
! h2 = occ(i,ispin)
|
||||
! if(list_orb_reverse_pert_rdm(h2).lt.0)return
|
||||
! h2 = list_orb_reverse_pert_rdm(h2)
|
||||
|
||||
! nkeys += 1
|
||||
! values(nkeys) = 0.5d0 * contrib * phase
|
||||
! keys(1,nkeys) = h1
|
||||
! keys(2,nkeys) = h2
|
||||
! keys(3,nkeys) = p1
|
||||
! keys(4,nkeys) = h2
|
||||
|
||||
! nkeys += 1
|
||||
! values(nkeys) = - 0.5d0 * contrib * phase
|
||||
! keys(1,nkeys) = h1
|
||||
! keys(2,nkeys) = h2
|
||||
! keys(3,nkeys) = h2
|
||||
! keys(4,nkeys) = p1
|
||||
!
|
||||
! nkeys += 1
|
||||
! values(nkeys) = 0.5d0 * contrib * phase
|
||||
! keys(1,nkeys) = h2
|
||||
! keys(2,nkeys) = h1
|
||||
! keys(3,nkeys) = h2
|
||||
! keys(4,nkeys) = p1
|
||||
|
||||
! nkeys += 1
|
||||
! values(nkeys) = - 0.5d0 * contrib * phase
|
||||
! keys(1,nkeys) = h2
|
||||
! keys(2,nkeys) = h1
|
||||
! keys(3,nkeys) = p1
|
||||
! keys(4,nkeys) = h2
|
||||
!enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
|
||||
implicit none
|
||||
integer, intent(in) :: sze_buff
|
||||
integer,intent(in) :: exc(0:2,2,2)
|
||||
double precision,intent(in) :: phase, contrib
|
||||
integer, intent(inout) :: nkeys, keys(4,sze_buff)
|
||||
double precision, intent(inout):: values(sze_buff)
|
||||
integer :: h1,h2,p1,p2
|
||||
|
||||
if (exc(0,1,1) == 1) then
|
||||
! Double alpha/beta
|
||||
h1 = exc(1,1,1)
|
||||
h2 = exc(1,1,2)
|
||||
p1 = exc(1,2,1)
|
||||
p2 = exc(1,2,2)
|
||||
! check if the orbitals involved are within the orbital range
|
||||
if(list_orb_reverse_pert_rdm(h1).lt.0)return
|
||||
h1 = list_orb_reverse_pert_rdm(h1)
|
||||
if(list_orb_reverse_pert_rdm(h2).lt.0)return
|
||||
h2 = list_orb_reverse_pert_rdm(h2)
|
||||
if(list_orb_reverse_pert_rdm(p1).lt.0)return
|
||||
p1 = list_orb_reverse_pert_rdm(p1)
|
||||
if(list_orb_reverse_pert_rdm(p2).lt.0)return
|
||||
p2 = list_orb_reverse_pert_rdm(p2)
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h1
|
||||
keys(2,nkeys) = h2
|
||||
keys(3,nkeys) = p1
|
||||
keys(4,nkeys) = p2
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = p1
|
||||
keys(2,nkeys) = p2
|
||||
keys(3,nkeys) = h1
|
||||
keys(4,nkeys) = h2
|
||||
|
||||
else
|
||||
if (exc(0,1,1) == 2) then
|
||||
! Double alpha/alpha
|
||||
h1 = exc(1,1,1)
|
||||
h2 = exc(2,1,1)
|
||||
p1 = exc(1,2,1)
|
||||
p2 = exc(2,2,1)
|
||||
else if (exc(0,1,2) == 2) then
|
||||
! Double beta
|
||||
h1 = exc(1,1,2)
|
||||
h2 = exc(2,1,2)
|
||||
p1 = exc(1,2,2)
|
||||
p2 = exc(2,2,2)
|
||||
endif
|
||||
! check if the orbitals involved are within the orbital range
|
||||
if(list_orb_reverse_pert_rdm(h1).lt.0)return
|
||||
h1 = list_orb_reverse_pert_rdm(h1)
|
||||
if(list_orb_reverse_pert_rdm(h2).lt.0)return
|
||||
h2 = list_orb_reverse_pert_rdm(h2)
|
||||
if(list_orb_reverse_pert_rdm(p1).lt.0)return
|
||||
p1 = list_orb_reverse_pert_rdm(p1)
|
||||
if(list_orb_reverse_pert_rdm(p2).lt.0)return
|
||||
p2 = list_orb_reverse_pert_rdm(p2)
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h1
|
||||
keys(2,nkeys) = h2
|
||||
keys(3,nkeys) = p1
|
||||
keys(4,nkeys) = p2
|
||||
|
||||
nkeys += 1
|
||||
values(nkeys) = - 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h1
|
||||
keys(2,nkeys) = h2
|
||||
keys(3,nkeys) = p2
|
||||
keys(4,nkeys) = p1
|
||||
|
||||
nkeys += 1
|
||||
values(nkeys) = 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h2
|
||||
keys(2,nkeys) = h1
|
||||
keys(3,nkeys) = p2
|
||||
keys(4,nkeys) = p1
|
||||
|
||||
nkeys += 1
|
||||
values(nkeys) = - 0.5d0 * contrib * phase
|
||||
keys(1,nkeys) = h2
|
||||
keys(2,nkeys) = h1
|
||||
keys(3,nkeys) = p1
|
||||
keys(4,nkeys) = p2
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
|
||||
@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data)
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
|
||||
PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
|
||||
PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max
|
||||
PROVIDE excitation_beta_max excitation_alpha_max excitation_max
|
||||
|
||||
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
|
||||
|
||||
|
||||
481
src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
Normal file
481
src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
Normal file
@ -0,0 +1,481 @@
|
||||
|
||||
subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
|
||||
use mmap_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR
|
||||
!
|
||||
! Dress_jj : DIAGONAL DRESSING of the Hamiltonian
|
||||
!
|
||||
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
|
||||
!
|
||||
! u_in : guess coefficients on the various states. Overwritten on exit
|
||||
!
|
||||
! sze : leftmost dimension of u_in
|
||||
!
|
||||
! sze : Number of determinants
|
||||
!
|
||||
! N_st : Number of eigenstates
|
||||
!
|
||||
! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze
|
||||
!
|
||||
! Initial guess vectors are not necessarily orthonormal
|
||||
!
|
||||
! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
|
||||
END_DOC
|
||||
integer, intent(in) :: sze, N_st, N_st_diag_in
|
||||
double precision, intent(in) :: H_jj(sze),Dress_jj(sze)
|
||||
double precision, intent(inout) :: u_in(sze,N_st_diag_in)
|
||||
double precision, intent(out) :: energies(N_st)
|
||||
external hcalc
|
||||
|
||||
integer :: iter, N_st_diag
|
||||
integer :: i,j,k,l,m
|
||||
logical, intent(inout) :: converged
|
||||
|
||||
double precision, external :: u_dot_v, u_dot_u
|
||||
|
||||
integer :: k_pairs, kl
|
||||
|
||||
integer :: iter2, itertot
|
||||
double precision, allocatable :: y(:,:), h(:,:), lambda(:)
|
||||
double precision, allocatable :: residual_norm(:)
|
||||
character*(16384) :: write_buffer
|
||||
double precision :: to_print(2,N_st)
|
||||
double precision :: cpu, wall
|
||||
integer :: shift, shift2, itermax, istate
|
||||
double precision :: r1, r2, alpha
|
||||
integer :: nproc_target
|
||||
integer :: order(N_st_diag_in)
|
||||
double precision :: cmax
|
||||
double precision, allocatable :: U(:,:), overlap(:,:)!, S_d(:,:)
|
||||
double precision, pointer :: W(:,:)
|
||||
logical :: disk_based
|
||||
double precision :: energy_shift(N_st_diag_in*davidson_sze_max)
|
||||
|
||||
include 'constants.include.F'
|
||||
|
||||
N_st_diag = N_st_diag_in
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
|
||||
if (N_st_diag*3 > sze) then
|
||||
print *, 'error in Davidson :'
|
||||
print *, 'Increase n_det_max_full to ', N_st_diag*3
|
||||
stop -1
|
||||
endif
|
||||
|
||||
itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
|
||||
itertot = 0
|
||||
|
||||
if (state_following) then
|
||||
allocate(overlap(N_st_diag*itermax, N_st_diag*itermax))
|
||||
else
|
||||
allocate(overlap(1,1)) ! avoid 'if' for deallocate
|
||||
endif
|
||||
overlap = 0.d0
|
||||
|
||||
provide threshold_davidson !nthreads_davidson
|
||||
call write_time(6)
|
||||
write(6,'(A)') ''
|
||||
write(6,'(A)') 'Davidson Diagonalization'
|
||||
write(6,'(A)') '------------------------'
|
||||
write(6,'(A)') ''
|
||||
|
||||
! Find max number of cores to fit in memory
|
||||
! -----------------------------------------
|
||||
|
||||
nproc_target = nproc
|
||||
double precision :: rss
|
||||
integer :: maxab
|
||||
maxab = sze
|
||||
|
||||
m=1
|
||||
disk_based = .False.
|
||||
call resident_memory(rss)
|
||||
do
|
||||
r1 = 8.d0 * &! bytes
|
||||
( dble(sze)*(N_st_diag*itermax) &! U
|
||||
+ 1.d0*dble(sze*m)*(N_st_diag*itermax) &! W
|
||||
+ 2.0d0*(N_st_diag*itermax)**2 &! h,y
|
||||
+ 2.d0*(N_st_diag*itermax) &! s2,lambda
|
||||
+ 1.d0*(N_st_diag) &! residual_norm
|
||||
! In H_S2_u_0_nstates_zmq
|
||||
+ 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on collector
|
||||
+ 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on slave
|
||||
+ 0.5d0*maxab &! idx0 in H_S2_u_0_nstates_openmp_work_*
|
||||
+ nproc_target * &! In OMP section
|
||||
( 1.d0*(N_int*maxab) &! buffer
|
||||
+ 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx
|
||||
) / 1024.d0**3
|
||||
|
||||
if (nproc_target == 0) then
|
||||
call check_mem(r1,irp_here)
|
||||
nproc_target = 1
|
||||
exit
|
||||
endif
|
||||
|
||||
if (r1+rss < qp_max_mem) then
|
||||
exit
|
||||
endif
|
||||
|
||||
if (itermax > 4) then
|
||||
itermax = itermax - 1
|
||||
else if (m==1.and.disk_based_davidson) then
|
||||
m=0
|
||||
disk_based = .True.
|
||||
itermax = 6
|
||||
else
|
||||
nproc_target = nproc_target - 1
|
||||
endif
|
||||
|
||||
enddo
|
||||
nthreads_davidson = nproc_target
|
||||
TOUCH nthreads_davidson
|
||||
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 basis functions')
|
||||
call write_int(6,nproc_target,'Number of threads for diagonalization')
|
||||
call write_double(6, r1, 'Memory(Gb)')
|
||||
if (disk_based) then
|
||||
print *, 'Using swap space to reduce RAM'
|
||||
endif
|
||||
|
||||
double precision, allocatable :: H_jj_tmp(:)
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (sze > 0)
|
||||
allocate(H_jj_tmp(sze))
|
||||
|
||||
do i=1,sze
|
||||
H_jj_tmp(i) = H_jj(i) + Dress_jj(i)
|
||||
enddo
|
||||
|
||||
!---------------
|
||||
|
||||
write(6,'(A)') ''
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = 'Iter'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' Energy Residual '
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
|
||||
|
||||
allocate(W(sze,N_st_diag*itermax))
|
||||
|
||||
allocate( &
|
||||
! Large
|
||||
U(sze,N_st_diag*itermax), &
|
||||
! Small
|
||||
h(N_st_diag*itermax,N_st_diag*itermax), &
|
||||
y(N_st_diag*itermax,N_st_diag*itermax), &
|
||||
residual_norm(N_st_diag), &
|
||||
lambda(N_st_diag*itermax))
|
||||
|
||||
h = 0.d0
|
||||
U = 0.d0
|
||||
y = 0.d0
|
||||
|
||||
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (N_st_diag >= N_st)
|
||||
ASSERT (sze > 0)
|
||||
|
||||
! Davidson iterations
|
||||
! ===================
|
||||
|
||||
converged = .False.
|
||||
|
||||
! Initialize from N_st to N_st_diat with gaussian random numbers
|
||||
! to be sure to have overlap with any eigenvectors
|
||||
do k=N_st+1,N_st_diag
|
||||
u_in(k,k) = 10.d0
|
||||
do i=1,sze
|
||||
call random_number(r1)
|
||||
call random_number(r2)
|
||||
r1 = dsqrt(-2.d0*dlog(r1))
|
||||
r2 = dtwo_pi*r2
|
||||
u_in(i,k) = r1*dcos(r2)
|
||||
enddo
|
||||
enddo
|
||||
! Normalize all states
|
||||
do k=1,N_st_diag
|
||||
call normalize(u_in(1,k),sze)
|
||||
enddo
|
||||
|
||||
! Copy from the guess input "u_in" to the working vectors "U"
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
do while (.not.converged)
|
||||
itertot = itertot+1
|
||||
if (itertot == 8) then
|
||||
exit
|
||||
endif
|
||||
|
||||
do iter=1,itermax-1
|
||||
|
||||
shift = N_st_diag*(iter-1)
|
||||
shift2 = N_st_diag*iter
|
||||
|
||||
if ((iter > 1).or.(itertot == 1)) then
|
||||
! Compute |W_k> = \sum_i |i><i|H|u_k>
|
||||
! -----------------------------------
|
||||
|
||||
! Gram-Schmidt to orthogonalize all new guess with the previous vectors
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag)
|
||||
! where sze is the size of the vector, N_st_diag is the number of states
|
||||
call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze)
|
||||
! Compute then the DIAGONAL PART OF THE DRESSING
|
||||
! <i|W_k> += Dress_jj(i) * <i|U>
|
||||
call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze)
|
||||
else
|
||||
! Already computed in update below
|
||||
continue
|
||||
endif
|
||||
|
||||
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
|
||||
! -------------------------------------------
|
||||
|
||||
call dgemm('T','N', shift2, shift2, sze, &
|
||||
1.d0, U, size(U,1), W, size(W,1), &
|
||||
0.d0, h, size(h,1))
|
||||
|
||||
! Diagonalize h y = lambda y
|
||||
! ---------------
|
||||
|
||||
call lapack_diag(lambda,y,h,size(h,1),shift2)
|
||||
|
||||
if (state_following) then
|
||||
|
||||
overlap = -1.d0
|
||||
do k=1,shift2
|
||||
do i=1,shift2
|
||||
overlap(k,i) = dabs(y(k,i))
|
||||
enddo
|
||||
enddo
|
||||
do k=1,N_st
|
||||
cmax = -1.d0
|
||||
do i=1,N_st
|
||||
if (overlap(i,k) > cmax) then
|
||||
cmax = overlap(i,k)
|
||||
order(k) = i
|
||||
endif
|
||||
enddo
|
||||
do i=1,N_st_diag
|
||||
overlap(order(k),i) = -1.d0
|
||||
enddo
|
||||
enddo
|
||||
overlap = y
|
||||
do k=1,N_st
|
||||
l = order(k)
|
||||
if (k /= l) then
|
||||
y(1:shift2,k) = overlap(1:shift2,l)
|
||||
endif
|
||||
enddo
|
||||
do k=1,N_st
|
||||
overlap(k,1) = lambda(k)
|
||||
enddo
|
||||
do k=1,N_st
|
||||
l = order(k)
|
||||
if (k /= l) then
|
||||
lambda(k) = overlap(l,1)
|
||||
endif
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! Express eigenvectors of h in the determinant basis
|
||||
! --------------------------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, &
|
||||
1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1))
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, &
|
||||
1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1))
|
||||
|
||||
! Compute residual vector and davidson step
|
||||
! -----------------------------------------
|
||||
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,shift2+k) = &
|
||||
(lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
|
||||
/max(H_jj_tmp(i) - lambda (k),1.d-2)
|
||||
enddo
|
||||
|
||||
if (k <= N_st) then
|
||||
residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
|
||||
to_print(1,k) = lambda(k)
|
||||
to_print(2,k) = residual_norm(k)
|
||||
endif
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
|
||||
if ((itertot>1).and.(iter == 1)) then
|
||||
!don't print
|
||||
continue
|
||||
else
|
||||
write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
|
||||
endif
|
||||
|
||||
! Check convergence
|
||||
if (iter > 1) then
|
||||
converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
|
||||
endif
|
||||
|
||||
|
||||
do k=1,N_st
|
||||
if (residual_norm(k) > 1.e8) then
|
||||
print *, 'Davidson failed'
|
||||
stop -1
|
||||
endif
|
||||
enddo
|
||||
if (converged) then
|
||||
exit
|
||||
endif
|
||||
|
||||
logical, external :: qp_stop
|
||||
if (qp_stop()) then
|
||||
converged = .True.
|
||||
exit
|
||||
endif
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
|
||||
W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
W(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
|
||||
U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
call ortho_qr(U,size(U,1),sze,N_st_diag)
|
||||
call ortho_qr(U,size(U,1),sze,N_st_diag)
|
||||
do j=1,N_st_diag
|
||||
k=1
|
||||
do while ((k<sze).and.(U(k,j) == 0.d0))
|
||||
k = k+1
|
||||
enddo
|
||||
if (U(k,j) * u_in(k,j) < 0.d0) then
|
||||
do i=1,sze
|
||||
W(i,j) = -W(i,j)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do k=1,N_st
|
||||
energies(k) = lambda(k)
|
||||
enddo
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') trim(write_buffer)
|
||||
write(6,'(A)') ''
|
||||
call write_time(6)
|
||||
|
||||
deallocate(W)
|
||||
|
||||
deallocate ( &
|
||||
residual_norm, &
|
||||
U, h, &
|
||||
y, &
|
||||
lambda &
|
||||
)
|
||||
deallocate(overlap)
|
||||
FREE nthreads_davidson
|
||||
end
|
||||
|
||||
subroutine hcalc_template(v,u,N_st,sze)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Template of routine for the application of H
|
||||
!
|
||||
! Here, it is done with the Hamiltonian matrix
|
||||
!
|
||||
! on the set of determinants of psi_det
|
||||
!
|
||||
! Computes $v = H | u \rangle$
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u(sze,N_st)
|
||||
double precision, intent(inout) :: v(sze,N_st)
|
||||
integer :: i,j,istate
|
||||
v = 0.d0
|
||||
do istate = 1, N_st
|
||||
do i = 1, sze
|
||||
do j = 1, sze
|
||||
v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate)
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, sze
|
||||
v(i,istate) += u(i,istate) * nuclear_repulsion
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Routine that computes the diagonal part of the dressing
|
||||
!
|
||||
! v(i) += u(i) * dress_diag(i)
|
||||
!
|
||||
! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized
|
||||
!
|
||||
! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE
|
||||
END_DOC
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u(sze,N_st),dress_diag(sze)
|
||||
double precision, intent(inout) :: v(sze,N_st)
|
||||
integer :: i,istate
|
||||
do istate = 1, N_st
|
||||
do i = 1, sze
|
||||
v(i,istate) += dress_diag(i) * u(i,istate)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
518
src/dav_general_mat/dav_double_dress_ext_rout.irp.f
Normal file
518
src/dav_general_mat/dav_double_dress_ext_rout.irp.f
Normal file
@ -0,0 +1,518 @@
|
||||
subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc)
|
||||
use mmap_module
|
||||
BEGIN_DOC
|
||||
! Generic Davidson diagonalization with TWO DRESSING VECTORS
|
||||
!
|
||||
! Dress_jj : DIAGONAL DRESSING of the Hamiltonian
|
||||
!
|
||||
! Dressing_vec : COLUMN / LINE DRESSING VECTOR
|
||||
!
|
||||
! idx_dress : position of the basis function used to use the Dressing_vec (usually the largest coeff)
|
||||
!
|
||||
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
|
||||
!
|
||||
! u_in : guess coefficients on the various states. Overwritten on exit
|
||||
!
|
||||
! sze : leftmost dimension of u_in
|
||||
!
|
||||
! sze : Number of determinants
|
||||
!
|
||||
! N_st : Number of eigenstates
|
||||
!
|
||||
! N_st_diag : Number of states in which H is diagonalized. Assumed > sze
|
||||
!
|
||||
! Initial guess vectors are not necessarily orthonormal
|
||||
!
|
||||
! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
|
||||
END_DOC
|
||||
implicit none
|
||||
integer, intent(in) :: sze, N_st, N_st_diag, idx_dress
|
||||
double precision, intent(in) :: H_jj(sze),Dress_jj(sze),Dressing_vec(sze,N_st)
|
||||
double precision, intent(inout) :: u_in(sze,N_st_diag)
|
||||
double precision, intent(out) :: energies(N_st_diag)
|
||||
logical, intent(out) :: converged
|
||||
external hcalc
|
||||
|
||||
double precision, allocatable :: H_jj_tmp(:)
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (sze > 0)
|
||||
allocate(H_jj_tmp(sze))
|
||||
|
||||
do i=1,sze
|
||||
H_jj_tmp(i) = H_jj(i) + Dress_jj(i)
|
||||
enddo
|
||||
do k=1,N_st
|
||||
do i=1,sze
|
||||
H_jj_tmp(i) += u_in(i,k) * Dressing_vec(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
integer :: iter
|
||||
integer :: i,j,k,l,m
|
||||
|
||||
double precision, external :: u_dot_v, u_dot_u
|
||||
|
||||
integer :: k_pairs, kl
|
||||
|
||||
integer :: iter2, itertot
|
||||
double precision, allocatable :: y(:,:), h(:,:), lambda(:)
|
||||
double precision, allocatable :: s_tmp(:,:)
|
||||
double precision, allocatable :: residual_norm(:),inv_c_idx_dress_vec(:)
|
||||
character*(16384) :: write_buffer
|
||||
double precision :: to_print(2,N_st),inv_c_idx_dress
|
||||
double precision :: cpu, wall
|
||||
integer :: shift, shift2, itermax, istate
|
||||
double precision :: r1, r2, alpha
|
||||
logical :: state_ok(N_st_diag*davidson_sze_max)
|
||||
integer :: nproc_target
|
||||
integer :: order(N_st_diag)
|
||||
double precision :: cmax
|
||||
double precision, allocatable :: U(:,:), overlap(:,:)
|
||||
double precision, pointer :: W(:,:)
|
||||
logical :: disk_based
|
||||
double precision :: energy_shift(N_st_diag*davidson_sze_max)
|
||||
|
||||
|
||||
allocate(inv_c_idx_dress_vec(N_st))
|
||||
inv_c_idx_dress = 1.d0/u_in(idx_dress,1)
|
||||
do i = 1, N_st
|
||||
inv_c_idx_dress_vec(i) = 1.d0/u_in(idx_dress,i)
|
||||
enddo
|
||||
include 'constants.include.F'
|
||||
|
||||
integer :: N_st_diag_in
|
||||
N_st_diag_in = N_st_diag
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
|
||||
if (N_st_diag_in*3 > sze) then
|
||||
print *, 'error in Davidson :'
|
||||
print *, 'Increase n_det_max_full to ', N_st_diag_in*3
|
||||
stop -1
|
||||
endif
|
||||
|
||||
itermax = max(2,min(davidson_sze_max, sze/N_st_diag_in))+1
|
||||
itertot = 0
|
||||
|
||||
if (state_following) then
|
||||
allocate(overlap(N_st_diag_in*itermax, N_st_diag_in*itermax))
|
||||
else
|
||||
allocate(overlap(1,1)) ! avoid 'if' for deallocate
|
||||
endif
|
||||
overlap = 0.d0
|
||||
|
||||
call write_time(6)
|
||||
write(6,'(A)') ''
|
||||
write(6,'(A)') 'Davidson Diagonalization'
|
||||
write(6,'(A)') '------------------------'
|
||||
write(6,'(A)') ''
|
||||
|
||||
! Find max number of cores to fit in memory
|
||||
! -----------------------------------------
|
||||
|
||||
nproc_target = nproc
|
||||
double precision :: rss
|
||||
integer :: maxab
|
||||
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
|
||||
|
||||
m=1
|
||||
disk_based = .False.
|
||||
call resident_memory(rss)
|
||||
do
|
||||
r1 = 8.d0 * &! bytes
|
||||
( dble(sze)*(N_st_diag_in*itermax) &! U
|
||||
+ 1.0d0*dble(sze*m)*(N_st_diag_in*itermax) &! W
|
||||
+ 3.0d0*(N_st_diag_in*itermax)**2 &! h,y,s_tmp
|
||||
+ 1.d0*(N_st_diag_in*itermax) &! lambda
|
||||
+ 1.d0*(N_st_diag_in) &! residual_norm
|
||||
! In H_u_0_nstates_zmq
|
||||
+ 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on collector
|
||||
+ 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on slave
|
||||
+ 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_*
|
||||
+ nproc_target * &! In OMP section
|
||||
( 1.d0*(N_int*maxab) &! buffer
|
||||
+ 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx
|
||||
) / 1024.d0**3
|
||||
|
||||
if (nproc_target == 0) then
|
||||
call check_mem(r1,irp_here)
|
||||
nproc_target = 1
|
||||
exit
|
||||
endif
|
||||
|
||||
if (r1+rss < qp_max_mem) then
|
||||
exit
|
||||
endif
|
||||
|
||||
if (itermax > 4) then
|
||||
itermax = itermax - 1
|
||||
else if (m==1.and.disk_based_davidson) then
|
||||
m=0
|
||||
disk_based = .True.
|
||||
itermax = 6
|
||||
else
|
||||
nproc_target = nproc_target - 1
|
||||
endif
|
||||
|
||||
enddo
|
||||
nthreads_davidson = nproc_target
|
||||
TOUCH nthreads_davidson
|
||||
call write_int(6,N_st,'Number of states')
|
||||
call write_int(6,N_st_diag_in,'Number of states in diagonalization')
|
||||
call write_int(6,sze,'Number of basis functions ')
|
||||
call write_int(6,nproc_target,'Number of threads for diagonalization')
|
||||
call write_double(6, r1, 'Memory(Gb)')
|
||||
if (disk_based) then
|
||||
print *, 'Using swap space to reduce RAM'
|
||||
endif
|
||||
|
||||
!---------------
|
||||
|
||||
write(6,'(A)') ''
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = 'Iter'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' Energy Residual '
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
|
||||
|
||||
allocate(W(sze,N_st_diag_in*itermax))
|
||||
|
||||
allocate( &
|
||||
! Large
|
||||
U(sze,N_st_diag_in*itermax), &
|
||||
|
||||
! Small
|
||||
h(N_st_diag_in*itermax,N_st_diag_in*itermax), &
|
||||
y(N_st_diag_in*itermax,N_st_diag_in*itermax), &
|
||||
s_tmp(N_st_diag_in*itermax,N_st_diag_in*itermax), &
|
||||
residual_norm(N_st_diag_in), &
|
||||
lambda(N_st_diag_in*itermax))
|
||||
|
||||
h = 0.d0
|
||||
U = 0.d0
|
||||
y = 0.d0
|
||||
s_tmp = 0.d0
|
||||
|
||||
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (N_st_diag_in >= N_st)
|
||||
ASSERT (sze > 0)
|
||||
|
||||
! Davidson iterations
|
||||
! ===================
|
||||
|
||||
converged = .False.
|
||||
|
||||
do k=N_st+1,N_st_diag_in
|
||||
do i=1,sze
|
||||
call random_number(r1)
|
||||
call random_number(r2)
|
||||
r1 = dsqrt(-2.d0*dlog(r1))
|
||||
r2 = dtwo_pi*r2
|
||||
u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
|
||||
enddo
|
||||
u_in(k,k) = u_in(k,k) + 10.d0
|
||||
enddo
|
||||
do k=1,N_st_diag_in
|
||||
call normalize(u_in(1,k),sze)
|
||||
enddo
|
||||
|
||||
do k=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
do while (.not.converged)
|
||||
itertot = itertot+1
|
||||
if (itertot == 2) then
|
||||
exit
|
||||
endif
|
||||
|
||||
do iter=1,itermax-1
|
||||
|
||||
shift = N_st_diag_in*(iter-1)
|
||||
shift2 = N_st_diag_in*iter
|
||||
|
||||
if ((iter > 1).or.(itertot == 1)) then
|
||||
! Compute |W_k> = \sum_i |i><i|H|u_k>
|
||||
! -----------------------------------
|
||||
call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag_in,sze)
|
||||
! Compute then the DIAGONAL PART OF THE DRESSING
|
||||
! <i|W_k> += Dress_jj(i) * <i|U>
|
||||
call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze)
|
||||
else
|
||||
! Already computed in update below
|
||||
continue
|
||||
endif
|
||||
|
||||
|
||||
if (N_st == 1) then
|
||||
|
||||
l = idx_dress
|
||||
double precision :: f
|
||||
f = inv_c_idx_dress
|
||||
do istate=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
W(i,shift+istate) += Dressing_vec(i,1) *f * U(l,shift+istate)
|
||||
W(l,shift+istate) += Dressing_vec(i,1) *f * U(i,shift+istate)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else
|
||||
print*,'dav_double_dressed routine not yet implemented for N_st > 1'
|
||||
!
|
||||
! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, &
|
||||
! psi_coef, size(psi_coef,1), &
|
||||
! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
|
||||
!
|
||||
! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, &
|
||||
! Dressing_vec, size(Dressing_vec,1), s_tmp, size(s_tmp,1), &
|
||||
! 1.d0, W(1,shift+1), size(W,1))
|
||||
!
|
||||
!
|
||||
! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, &
|
||||
! Dressing_vec, size(Dressing_vec,1), &
|
||||
! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
|
||||
!
|
||||
! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, &
|
||||
! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
|
||||
! 1.d0, W(1,shift+1), size(W,1))
|
||||
!
|
||||
endif
|
||||
|
||||
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
|
||||
! -------------------------------------------
|
||||
|
||||
call dgemm('T','N', shift2, shift2, sze, &
|
||||
1.d0, U, size(U,1), W, size(W,1), &
|
||||
0.d0, h, size(h,1))
|
||||
call dgemm('T','N', shift2, shift2, sze, &
|
||||
1.d0, U, size(U,1), U, size(U,1), &
|
||||
0.d0, s_tmp, size(s_tmp,1))
|
||||
|
||||
! Diagonalize h
|
||||
! ---------------
|
||||
|
||||
integer :: lwork, info
|
||||
double precision, allocatable :: work(:)
|
||||
|
||||
y = h
|
||||
lwork = -1
|
||||
allocate(work(1))
|
||||
call dsygv(1,'V','U',shift2,y,size(y,1), &
|
||||
s_tmp,size(s_tmp,1), lambda, work,lwork,info)
|
||||
lwork = int(work(1))
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
call dsygv(1,'V','U',shift2,y,size(y,1), &
|
||||
s_tmp,size(s_tmp,1), lambda, work,lwork,info)
|
||||
deallocate(work)
|
||||
if (info /= 0) then
|
||||
stop 'DSYGV Diagonalization failed'
|
||||
endif
|
||||
|
||||
! Compute Energy for each eigenvector
|
||||
! -----------------------------------
|
||||
|
||||
call dgemm('N','N',shift2,shift2,shift2, &
|
||||
1.d0, h, size(h,1), y, size(y,1), &
|
||||
0.d0, s_tmp, size(s_tmp,1))
|
||||
|
||||
call dgemm('T','N',shift2,shift2,shift2, &
|
||||
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
|
||||
0.d0, h, size(h,1))
|
||||
|
||||
do k=1,shift2
|
||||
lambda(k) = h(k,k)
|
||||
enddo
|
||||
|
||||
if (state_following) then
|
||||
|
||||
overlap = -1.d0
|
||||
do k=1,shift2
|
||||
do i=1,shift2
|
||||
overlap(k,i) = dabs(y(k,i))
|
||||
enddo
|
||||
enddo
|
||||
do k=1,N_st
|
||||
cmax = -1.d0
|
||||
do i=1,N_st
|
||||
if (overlap(i,k) > cmax) then
|
||||
cmax = overlap(i,k)
|
||||
order(k) = i
|
||||
endif
|
||||
enddo
|
||||
do i=1,N_st_diag_in
|
||||
overlap(order(k),i) = -1.d0
|
||||
enddo
|
||||
enddo
|
||||
overlap = y
|
||||
do k=1,N_st
|
||||
l = order(k)
|
||||
if (k /= l) then
|
||||
y(1:shift2,k) = overlap(1:shift2,l)
|
||||
endif
|
||||
enddo
|
||||
do k=1,N_st
|
||||
overlap(k,1) = lambda(k)
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! Express eigenvectors of h in the determinant basis
|
||||
! --------------------------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag_in, shift2, &
|
||||
1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1))
|
||||
call dgemm('N','N', sze, N_st_diag_in, shift2, &
|
||||
1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1))
|
||||
|
||||
! Compute residual vector and davidson step
|
||||
! -----------------------------------------
|
||||
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
|
||||
do k=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
U(i,shift2+k) = &
|
||||
(lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
|
||||
/max(H_jj_tmp(i) - lambda (k),1.d-2)
|
||||
enddo
|
||||
|
||||
if (k <= N_st) then
|
||||
residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
|
||||
to_print(1,k) = lambda(k)
|
||||
to_print(2,k) = residual_norm(k)
|
||||
endif
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
|
||||
if ((itertot>1).and.(iter == 1)) then
|
||||
!don't print
|
||||
continue
|
||||
else
|
||||
write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
|
||||
endif
|
||||
|
||||
! Check convergence
|
||||
if (iter > 1) then
|
||||
converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
|
||||
endif
|
||||
|
||||
do k=1,N_st
|
||||
if (residual_norm(k) > 1.d8) then
|
||||
print *, 'Davidson failed'
|
||||
stop -1
|
||||
endif
|
||||
enddo
|
||||
if (converged) then
|
||||
exit
|
||||
endif
|
||||
|
||||
logical, external :: qp_stop
|
||||
if (qp_stop()) then
|
||||
converged = .True.
|
||||
exit
|
||||
endif
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
! Re-contract U and update W
|
||||
! --------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, &
|
||||
W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
do k=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
W(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, &
|
||||
U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
|
||||
do k=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
call nullify_small_elements(sze,N_st_diag_in,U,size(U,1),threshold_davidson_pt2)
|
||||
do k=1,N_st_diag_in
|
||||
do i=1,sze
|
||||
u_in(i,k) = U(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do k=1,N_st_diag_in
|
||||
energies(k) = lambda(k)
|
||||
enddo
|
||||
write_buffer = '======'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') trim(write_buffer)
|
||||
write(6,'(A)') ''
|
||||
call write_time(6)
|
||||
|
||||
deallocate(W)
|
||||
|
||||
deallocate ( &
|
||||
residual_norm, &
|
||||
U, overlap, &
|
||||
h, y, s_tmp, &
|
||||
lambda &
|
||||
)
|
||||
FREE nthreads_davidson
|
||||
end
|
||||
|
||||
|
||||
subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Routine that computes the diagonal part of the dressing
|
||||
!
|
||||
! v(i) += u(i) * dress_diag(i)
|
||||
!
|
||||
! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized
|
||||
!
|
||||
! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE
|
||||
END_DOC
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u(sze,N_st),dress_diag(sze)
|
||||
double precision, intent(inout) :: v(sze,N_st)
|
||||
integer :: i,istate
|
||||
do istate = 1, N_st
|
||||
do i = 1, sze
|
||||
v(i,istate) += dress_diag(i) * u(i,istate)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
485
src/dav_general_mat/dav_dressed_ext_rout.irp.f
Normal file
485
src/dav_general_mat/dav_dressed_ext_rout.irp.f
Normal file
@ -0,0 +1,485 @@
|
||||
subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag,dressing_state,dressing_vec,idress,converged,hcalc)
|
||||
use mmap_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Davidson diagonalization.
|
||||
!
|
||||
! u_in : guess coefficients on the various states. Overwritten
|
||||
! on exit
|
||||
!
|
||||
! sze : leftmost dimension of u_in
|
||||
!
|
||||
! sze : Number of determinants
|
||||
!
|
||||
! N_st : Number of eigenstates
|
||||
!
|
||||
! Initial guess vectors are not necessarily orthonormal
|
||||
END_DOC
|
||||
integer, intent(in) :: sze, N_st, N_st_diag,idress
|
||||
double precision, intent(inout) :: u_in(sze,N_st_diag)
|
||||
double precision, intent(inout) :: H_jj(sze)
|
||||
double precision, intent(out) :: energies(N_st_diag)
|
||||
double precision, intent(in) :: dressing_vec(sze,N_st)
|
||||
integer, intent(in) :: dressing_state
|
||||
logical, intent(out) :: converged
|
||||
external hcalc
|
||||
|
||||
double precision :: f
|
||||
|
||||
integer :: iter
|
||||
integer :: i,j,k,l,m
|
||||
|
||||
double precision, external :: u_dot_v, u_dot_u
|
||||
|
||||
integer :: k_pairs, kl
|
||||
|
||||
integer :: iter2, itertot
|
||||
double precision, allocatable :: y(:,:), h(:,:), lambda(:)
|
||||
double precision, allocatable :: s_tmp(:,:)
|
||||
double precision :: diag_h_mat_elem
|
||||
double precision, allocatable :: residual_norm(:)
|
||||
character*(16384) :: write_buffer
|
||||
double precision :: to_print(2,N_st)
|
||||
double precision :: cpu, wall
|
||||
integer :: shift, shift2, itermax, istate
|
||||
double precision :: r1, r2, alpha
|
||||
logical :: state_ok(N_st_diag*davidson_sze_max)
|
||||
integer :: nproc_target
|
||||
integer :: order(N_st_diag)
|
||||
double precision :: cmax
|
||||
double precision, allocatable :: U(:,:), overlap(:,:)
|
||||
double precision, pointer :: W(:,:)
|
||||
logical :: disk_based
|
||||
double precision :: energy_shift(N_st_diag*davidson_sze_max)
|
||||
!!!! TO CHANGE !!!!
|
||||
integer :: idx_dress(1)
|
||||
idx_dress = idress
|
||||
|
||||
|
||||
if (dressing_state > 0) then
|
||||
do k=1,N_st
|
||||
do i=1,sze
|
||||
H_jj(i) += u_in(i,k) * dressing_vec(i,k)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
l = idx_dress(1)
|
||||
f = 1.0d0/u_in(l,1)
|
||||
|
||||
include 'constants.include.F'
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
|
||||
if (N_st_diag*3 > sze) then
|
||||
print *, 'error in Davidson :'
|
||||
print *, 'Increase n_det_max_full to ', N_st_diag*3
|
||||
stop -1
|
||||
endif
|
||||
|
||||
itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
|
||||
itertot = 0
|
||||
|
||||
if (state_following) then
|
||||
allocate(overlap(N_st_diag*itermax, N_st_diag*itermax))
|
||||
else
|
||||
allocate(overlap(1,1)) ! avoid 'if' for deallocate
|
||||
endif
|
||||
overlap = 0.d0
|
||||
|
||||
|
||||
provide threshold_davidson !nthreads_davidson
|
||||
call write_time(6)
|
||||
write(6,'(A)') ''
|
||||
write(6,'(A)') 'Davidson Diagonalization'
|
||||
write(6,'(A)') '------------------------'
|
||||
write(6,'(A)') ''
|
||||
|
||||
! Find max number of cores to fit in memory
|
||||
! -----------------------------------------
|
||||
|
||||
nproc_target = nproc
|
||||
double precision :: rss
|
||||
integer :: maxab
|
||||
! maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
|
||||
maxab = sze
|
||||
|
||||
m=1
|
||||
disk_based = .False.
|
||||
call resident_memory(rss)
|
||||
do
|
||||
r1 = 8.d0 * &! bytes
|
||||
( dble(sze)*(N_st_diag*itermax) &! U
|
||||
+ 1.0d0*dble(sze*m)*(N_st_diag*itermax) &! W
|
||||
+ 3.0d0*(N_st_diag*itermax)**2 &! h,y,s_tmp
|
||||
+ 1.d0*(N_st_diag*itermax) &! lambda
|
||||
+ 1.d0*(N_st_diag) &! residual_norm
|
||||
! In H_u_0_nstates_zmq
|
||||
+ 2.d0*(N_st_diag*N_det) &! u_t, v_t, on collector
|
||||
+ 2.d0*(N_st_diag*N_det) &! u_t, v_t, on slave
|
||||
+ 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_*
|
||||
+ nproc_target * &! In OMP section
|
||||
( 1.d0*(N_int*maxab) &! buffer
|
||||
+ 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx
|
||||
) / 1024.d0**3
|
||||
|
||||
if (nproc_target == 0) then
|
||||
call check_mem(r1,irp_here)
|
||||
nproc_target = 1
|
||||
exit
|
||||
endif
|
||||
|
||||
if (r1+rss < qp_max_mem) then
|
||||
exit
|
||||
endif
|
||||
|
||||
if (itermax > 4) then
|
||||
itermax = itermax - 1
|
||||
else if (m==1.and.disk_based_davidson) then
|
||||
m=0
|
||||
disk_based = .True.
|
||||
itermax = 6
|
||||
else
|
||||
nproc_target = nproc_target - 1
|
||||
endif
|
||||
|
||||
enddo
|
||||
nthreads_davidson = nproc_target
|
||||
TOUCH nthreads_davidson
|
||||
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 basis function')
|
||||
call write_int(6,nproc_target,'Number of threads for diagonalization')
|
||||
call write_double(6, r1, 'Memory(Gb)')
|
||||
if (disk_based) then
|
||||
print *, 'Using swap space to reduce RAM'
|
||||
endif
|
||||
|
||||
!---------------
|
||||
|
||||
write(6,'(A)') ''
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = 'Iter'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' Energy Residual '
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
write_buffer = '====='
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
|
||||
|
||||
allocate(W(sze,N_st_diag*itermax))
|
||||
|
||||
allocate( &
|
||||
! Large
|
||||
U(sze,N_st_diag*itermax), &
|
||||
|
||||
! Small
|
||||
h(N_st_diag*itermax,N_st_diag*itermax), &
|
||||
y(N_st_diag*itermax,N_st_diag*itermax), &
|
||||
s_tmp(N_st_diag*itermax,N_st_diag*itermax), &
|
||||
residual_norm(N_st_diag), &
|
||||
lambda(N_st_diag*itermax))
|
||||
|
||||
h = 0.d0
|
||||
U = 0.d0
|
||||
y = 0.d0
|
||||
s_tmp = 0.d0
|
||||
|
||||
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (N_st_diag >= N_st)
|
||||
ASSERT (sze > 0)
|
||||
|
||||
! Davidson iterations
|
||||
! ===================
|
||||
|
||||
converged = .False.
|
||||
|
||||
do k=N_st+1,N_st_diag
|
||||
do i=1,sze
|
||||
call random_number(r1)
|
||||
call random_number(r2)
|
||||
r1 = dsqrt(-2.d0*dlog(r1))
|
||||
r2 = dtwo_pi*r2
|
||||
u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
|
||||
enddo
|
||||
u_in(k,k) = u_in(k,k) + 10.d0
|
||||
enddo
|
||||
! Normalize all states
|
||||
do k=1,N_st_diag
|
||||
call normalize(u_in(1,k),sze)
|
||||
enddo
|
||||
! Copy from the guess input "u_in" to the working vectors "U"
|
||||
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
do while (.not.converged)
|
||||
itertot = itertot+1
|
||||
if (itertot == 8) then
|
||||
exit
|
||||
endif
|
||||
|
||||
do iter=1,itermax-1
|
||||
|
||||
shift = N_st_diag*(iter-1)
|
||||
shift2 = N_st_diag*iter
|
||||
|
||||
if ((iter > 1).or.(itertot == 1)) then
|
||||
! Compute |W_k> = \sum_i |i><i|H|u_k>
|
||||
! -----------------------------------
|
||||
! Gram-Schmidt to orthogonalize all new guess with the previous vectors
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag)
|
||||
! where sze is the size of the vector, N_st_diag is the number of states
|
||||
call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze)
|
||||
else
|
||||
! Already computed in update below
|
||||
continue
|
||||
endif
|
||||
|
||||
if (dressing_state > 0) then
|
||||
|
||||
if (N_st == 1) then
|
||||
|
||||
|
||||
do istate=1,N_st_diag
|
||||
do i=1,sze
|
||||
W(i,shift+istate) += dressing_vec(i,1) *f * U(l,shift+istate)
|
||||
W(l,shift+istate) += dressing_vec(i,1) *f * U(i,shift+istate)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
else
|
||||
print*,'Not implemented yet for multi state ...'
|
||||
stop
|
||||
! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
|
||||
! psi_coef, size(psi_coef,1), &
|
||||
! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
|
||||
!
|
||||
! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, &
|
||||
! dressing_vec, size(dressing_vec,1), s_tmp, size(s_tmp,1), &
|
||||
! 1.d0, W(1,shift+1), size(W,1))
|
||||
!
|
||||
!
|
||||
! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
|
||||
! dressing_vec, size(dressing_vec,1), &
|
||||
! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
|
||||
!
|
||||
! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, &
|
||||
! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
|
||||
! 1.d0, W(1,shift+1), size(W,1))
|
||||
|
||||
endif
|
||||
endif
|
||||
|
||||
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
|
||||
! -------------------------------------------
|
||||
|
||||
call dgemm('T','N', shift2, shift2, sze, &
|
||||
1.d0, U, size(U,1), W, size(W,1), &
|
||||
0.d0, h, size(h,1))
|
||||
call dgemm('T','N', shift2, shift2, sze, &
|
||||
1.d0, U, size(U,1), U, size(U,1), &
|
||||
0.d0, s_tmp, size(s_tmp,1))
|
||||
|
||||
! Diagonalize h
|
||||
! ---------------
|
||||
|
||||
integer :: lwork, info
|
||||
double precision, allocatable :: work(:)
|
||||
|
||||
y = h
|
||||
lwork = -1
|
||||
allocate(work(1))
|
||||
call dsygv(1,'V','U',shift2,y,size(y,1), &
|
||||
s_tmp,size(s_tmp,1), lambda, work,lwork,info)
|
||||
lwork = int(work(1))
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
call dsygv(1,'V','U',shift2,y,size(y,1), &
|
||||
s_tmp,size(s_tmp,1), lambda, work,lwork,info)
|
||||
deallocate(work)
|
||||
if (info /= 0) then
|
||||
stop 'DSYGV Diagonalization failed'
|
||||
endif
|
||||
|
||||
! Compute Energy for each eigenvector
|
||||
! -----------------------------------
|
||||
|
||||
call dgemm('N','N',shift2,shift2,shift2, &
|
||||
1.d0, h, size(h,1), y, size(y,1), &
|
||||
0.d0, s_tmp, size(s_tmp,1))
|
||||
|
||||
call dgemm('T','N',shift2,shift2,shift2, &
|
||||
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
|
||||
0.d0, h, size(h,1))
|
||||
|
||||
do k=1,shift2
|
||||
lambda(k) = h(k,k)
|
||||
enddo
|
||||
|
||||
if (state_following) then
|
||||
|
||||
overlap = -1.d0
|
||||
do k=1,shift2
|
||||
do i=1,shift2
|
||||
overlap(k,i) = dabs(y(k,i))
|
||||
enddo
|
||||
enddo
|
||||
do k=1,N_st
|
||||
cmax = -1.d0
|
||||
do i=1,N_st
|
||||
if (overlap(i,k) > cmax) then
|
||||
cmax = overlap(i,k)
|
||||
order(k) = i
|
||||
endif
|
||||
enddo
|
||||
do i=1,N_st_diag
|
||||
overlap(order(k),i) = -1.d0
|
||||
enddo
|
||||
enddo
|
||||
overlap = y
|
||||
do k=1,N_st
|
||||
l = order(k)
|
||||
if (k /= l) then
|
||||
y(1:shift2,k) = overlap(1:shift2,l)
|
||||
endif
|
||||
enddo
|
||||
do k=1,N_st
|
||||
overlap(k,1) = lambda(k)
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! Express eigenvectors of h in the determinant basis
|
||||
! --------------------------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, &
|
||||
1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1))
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, &
|
||||
1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1))
|
||||
|
||||
! Compute residual vector and davidson step
|
||||
! -----------------------------------------
|
||||
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,shift2+k) = &
|
||||
(lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
|
||||
/max(H_jj(i) - lambda (k),1.d-2)
|
||||
enddo
|
||||
|
||||
if (k <= N_st) then
|
||||
residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
|
||||
to_print(1,k) = lambda(k)
|
||||
to_print(2,k) = residual_norm(k)
|
||||
endif
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
|
||||
if ((itertot>1).and.(iter == 1)) then
|
||||
!don't print
|
||||
continue
|
||||
else
|
||||
write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
|
||||
endif
|
||||
|
||||
! Check convergence
|
||||
if (iter > 1) then
|
||||
converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
|
||||
endif
|
||||
|
||||
do k=1,N_st
|
||||
if (residual_norm(k) > 1.d8) then
|
||||
print *, 'Davidson failed'
|
||||
stop -1
|
||||
endif
|
||||
enddo
|
||||
if (converged) then
|
||||
exit
|
||||
endif
|
||||
|
||||
logical, external :: qp_stop
|
||||
if (qp_stop()) then
|
||||
converged = .True.
|
||||
exit
|
||||
endif
|
||||
|
||||
|
||||
enddo
|
||||
|
||||
! Re-contract U and update W
|
||||
! --------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
|
||||
W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
W(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
|
||||
U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
U(i,k) = u_in(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
call nullify_small_elements(sze,N_st_diag,U,size(U,1),threshold_davidson_pt2)
|
||||
do k=1,N_st_diag
|
||||
do i=1,sze
|
||||
u_in(i,k) = U(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do k=1,N_st_diag
|
||||
energies(k) = lambda(k)
|
||||
enddo
|
||||
write_buffer = '======'
|
||||
do i=1,N_st
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') trim(write_buffer)
|
||||
write(6,'(A)') ''
|
||||
call write_time(6)
|
||||
|
||||
deallocate(W)
|
||||
|
||||
deallocate ( &
|
||||
residual_norm, &
|
||||
U, overlap, &
|
||||
h, y, s_tmp, &
|
||||
lambda &
|
||||
)
|
||||
FREE nthreads_davidson
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,15 +1,15 @@
|
||||
|
||||
subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc)
|
||||
subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
|
||||
use mmap_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Davidson diagonalization with specific diagonal elements of the H matrix
|
||||
! Generic Davidson diagonalization
|
||||
!
|
||||
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
|
||||
!
|
||||
! u_in : guess coefficients on the various states. Overwritten on exit
|
||||
!
|
||||
! dim_in : leftmost dimension of u_in
|
||||
! sze : leftmost dimension of u_in
|
||||
!
|
||||
! sze : Number of determinants
|
||||
!
|
||||
@ -21,9 +21,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
||||
!
|
||||
! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
|
||||
END_DOC
|
||||
integer, intent(in) :: dim_in, sze, N_st, N_st_diag_in
|
||||
integer, intent(in) :: sze, N_st, N_st_diag_in
|
||||
double precision, intent(in) :: H_jj(sze)
|
||||
double precision, intent(inout) :: u_in(dim_in,N_st_diag_in)
|
||||
double precision, intent(inout) :: u_in(sze,N_st_diag_in)
|
||||
double precision, intent(out) :: energies(N_st)
|
||||
external hcalc
|
||||
|
||||
@ -157,19 +157,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
||||
write(6,'(A)') write_buffer(1:6+41*N_st)
|
||||
|
||||
|
||||
! if (disk_based) then
|
||||
! ! Create memory-mapped files for W and S
|
||||
! type(c_ptr) :: ptr_w, ptr_s
|
||||
! integer :: fd_s, fd_w
|
||||
! call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
|
||||
! 8, fd_w, .False., ptr_w)
|
||||
! call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),&
|
||||
! 4, fd_s, .False., ptr_s)
|
||||
! call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
|
||||
! call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/))
|
||||
! else
|
||||
allocate(W(sze,N_st_diag*itermax))
|
||||
! endif
|
||||
|
||||
allocate( &
|
||||
! Large
|
||||
@ -233,7 +221,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
||||
if ((iter > 1).or.(itertot == 1)) then
|
||||
! Compute |W_k> = \sum_i |i><i|H|u_k>
|
||||
! -----------------------------------
|
||||
|
||||
! Gram-Schmidt to orthogonalize all new guess with the previous vectors
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
call ortho_qr(U,size(U,1),sze,shift2)
|
||||
@ -357,6 +344,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
|
||||
|
||||
enddo
|
||||
|
||||
! Re-contract U and update W
|
||||
! --------------------------------
|
||||
|
||||
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
|
||||
W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
|
||||
do k=1,N_st_diag
|
||||
|
||||
@ -8,12 +8,13 @@ program test_dav
|
||||
touch read_wf
|
||||
PROVIDE threshold_davidson nthreads_davidson
|
||||
call routine
|
||||
call test_dav_dress
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
double precision, allocatable :: u_in(:,:), H_jj(:), energies(:),h_mat(:,:)
|
||||
integer :: dim_in,sze,N_st,N_st_diag_in,dressing_state
|
||||
integer :: dim_in,sze,N_st,N_st_diag_in
|
||||
logical :: converged
|
||||
integer :: i,j
|
||||
external hcalc_template
|
||||
@ -21,9 +22,8 @@ subroutine routine
|
||||
N_st_diag_in = N_states_diag
|
||||
sze = N_det
|
||||
dim_in = sze
|
||||
dressing_state = 0
|
||||
!!!! MARK THAT u_in mut dimensioned with "N_st_diag_in" as a second dimension
|
||||
allocate(u_in(dim_in,N_st_diag_in),H_jj(sze),h_mat(sze,sze),energies(N_st))
|
||||
allocate(u_in(dim_in,N_st_diag_in),H_jj(sze),h_mat(sze,sze),energies(N_st_diag_in))
|
||||
u_in = 0.d0
|
||||
do i = 1, N_st
|
||||
u_in(1,i) = 1.d0
|
||||
@ -42,7 +42,38 @@ subroutine routine
|
||||
print*,'energies = ',energies
|
||||
!!! hcalc_template is the routine that computes v = H u
|
||||
!!! and you can use the routine "davidson_general_ext_rout"
|
||||
call davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc_template)
|
||||
call davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc_template)
|
||||
print*,'energies = ',energies
|
||||
end
|
||||
|
||||
|
||||
subroutine test_dav_dress
|
||||
implicit none
|
||||
double precision, allocatable :: u_in(:,:), H_jj(:), energies(:)
|
||||
integer :: sze,N_st,N_st_diag_in,dressing_state
|
||||
logical :: converged
|
||||
integer :: i,j
|
||||
external hcalc_template
|
||||
double precision, allocatable :: dressing_vec(:)
|
||||
integer :: idress
|
||||
N_st = N_states
|
||||
N_st_diag_in = N_states_diag
|
||||
sze = N_det
|
||||
dressing_state = 0
|
||||
idress = 1
|
||||
!!!! MARK THAT u_in mut dimensioned with "N_st_diag_in" as a second dimension
|
||||
allocate(u_in(sze,N_st_diag_in),H_jj(sze),energies(N_st_diag_in))
|
||||
allocate(dressing_vec(sze))
|
||||
dressing_vec = 0.d0
|
||||
u_in = 0.d0
|
||||
do i = 1, N_st
|
||||
u_in(1,i) = 1.d0
|
||||
enddo
|
||||
do i = 1, sze
|
||||
H_jj(i) = H_matrix_all_dets(i,i) + nuclear_repulsion
|
||||
enddo
|
||||
print*,'dressing davidson '
|
||||
call davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag_in,dressing_state,dressing_vec,idress,converged,hcalc_template)
|
||||
print*,'energies(1) = ',energies(1)
|
||||
|
||||
end
|
||||
|
||||
@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
|
||||
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_energy_dressed(j),'Energy dressed of state '//trim(st))
|
||||
call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
|
||||
enddo
|
||||
|
||||
|
||||
@ -140,5 +140,4 @@ default: False
|
||||
type: Threshold
|
||||
doc: Cut-off to apply to the CI coefficients when the wave function is stored
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.
|
||||
|
||||
default: 1.e-14
|
||||
|
||||
@ -268,6 +268,44 @@ subroutine set_natural_mos
|
||||
soft_touch mo_occ
|
||||
end
|
||||
|
||||
|
||||
subroutine save_natural_mos_canon_label
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Save natural orbitals, obtained by diagonalization of the one-body density matrix in
|
||||
! the |MO| basis
|
||||
END_DOC
|
||||
call set_natural_mos_canon_label
|
||||
call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10)
|
||||
call orthonormalize_mos
|
||||
call save_mos
|
||||
end
|
||||
|
||||
subroutine set_natural_mos_canon_label
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Set natural orbitals, obtained by diagonalization of the one-body density matrix
|
||||
! in the |MO| basis
|
||||
END_DOC
|
||||
character*(64) :: label
|
||||
double precision, allocatable :: tmp(:,:)
|
||||
|
||||
label = "Canonical"
|
||||
integer :: i,j,iorb,jorb
|
||||
do i = 1, n_virt_orb
|
||||
iorb = list_virt(i)
|
||||
do j = 1, n_core_inact_act_orb
|
||||
jorb = list_core_inact_act(j)
|
||||
enddo
|
||||
enddo
|
||||
call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label)
|
||||
soft_touch mo_occ
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine set_natorb_no_ov_rot
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
||||
@ -438,7 +438,7 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gives the inidices(+1) of the bits set to 1 in the bit string
|
||||
! Gives the indices(+1) of the bits set to 1 in the bit string
|
||||
! For alpha/beta determinants.
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
@ -472,6 +472,35 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
|
||||
|
||||
end
|
||||
|
||||
!subroutine bitstring_to_list( string, list, n_elements, Nint)
|
||||
! use bitmasks
|
||||
! implicit none
|
||||
! BEGIN_DOC
|
||||
! ! Gives the indices(+1) of the bits set to 1 in the bit string
|
||||
! END_DOC
|
||||
! integer, intent(in) :: Nint
|
||||
! integer(bit_kind), intent(in) :: string(Nint)
|
||||
! integer, intent(out) :: list(Nint*bit_kind_size)
|
||||
! integer, intent(out) :: n_elements
|
||||
!
|
||||
! integer :: i, j, ishift
|
||||
! integer(bit_kind) :: l
|
||||
!
|
||||
! n_elements = 0
|
||||
! ishift = 1
|
||||
! do i=1,Nint
|
||||
! l = string(i)
|
||||
! do while (l /= 0_bit_kind)
|
||||
! j = trailz(l)
|
||||
! n_elements = n_elements + 1
|
||||
! l = ibclr(l,j)
|
||||
! list(n_elements) = ishift+j
|
||||
! enddo
|
||||
! ishift = ishift + bit_kind_size
|
||||
! enddo
|
||||
!
|
||||
!end
|
||||
|
||||
|
||||
subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
|
||||
use bitmasks
|
||||
|
||||
@ -8,6 +8,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
|
||||
double precision :: hij
|
||||
integer :: degree(N_det),idx(0:N_det)
|
||||
call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
|
||||
print*,'Providing the H_matrix_all_dets ...'
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) &
|
||||
!$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets)
|
||||
do i =1,N_det
|
||||
@ -18,6 +19,26 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
print*,'H_matrix_all_dets done '
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! |H| matrix on the basis of the Slater determinants defined by psi_det
|
||||
END_DOC
|
||||
integer :: i
|
||||
double precision :: hij
|
||||
integer :: degree(N_det)
|
||||
call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) &
|
||||
!$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets)
|
||||
do i =1,N_det
|
||||
call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
|
||||
H_matrix_diag_all_dets(i) = hij
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
@ -16,3 +16,8 @@ doc: Percentage of HF exchange in the DFT model
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.
|
||||
|
||||
[mu_dft_type]
|
||||
type: character*(32)
|
||||
doc: type of mu(r) for rsdft [ cst ]
|
||||
interface: ezfio, provider, ocaml
|
||||
default: cst
|
||||
|
||||
@ -6,3 +6,4 @@ ao_one_e_ints
|
||||
ao_two_e_ints
|
||||
mo_two_e_erf_ints
|
||||
ao_two_e_erf_ints
|
||||
mu_of_r
|
||||
|
||||
@ -8,3 +8,73 @@ BEGIN_PROVIDER [double precision, mu_erf_dft]
|
||||
mu_erf_dft = mu_erf
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)]
|
||||
implicit none
|
||||
integer :: i
|
||||
if(mu_dft_type == "Read")then
|
||||
call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_dft)
|
||||
else
|
||||
do i = 1, n_points_final_grid
|
||||
if(mu_dft_type == "cst")then
|
||||
mu_of_r_dft(i) = mu_erf_dft
|
||||
else if(mu_dft_type == "hf")then
|
||||
mu_of_r_dft(i) = mu_of_r_hf(i)
|
||||
else if(mu_dft_type == "rsc")then
|
||||
mu_of_r_dft(i) = mu_rsc_of_r(i)
|
||||
else if(mu_dft_type == "grad_rho")then
|
||||
mu_of_r_dft(i) = mu_grad_rho(i)
|
||||
else
|
||||
print*,'mu_dft_type is not of good type = ',mu_dft_type
|
||||
print*,'it must be of type Read, cst, hf, rsc'
|
||||
print*,'Stopping ...'
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)]
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: mu_rs_c,rho,r(3), dm_a, dm_b
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
|
||||
rho = dm_a + dm_b
|
||||
mu_rsc_of_r(i) = mu_rs_c(rho)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, mu_grad_rho, (n_points_final_grid)]
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: mu_grad_rho_func, r(3)
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
mu_grad_rho(i) = mu_grad_rho_func(r)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, mu_of_r_dft_average]
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: mu_rs_c,rho,r(3), dm_a, dm_b
|
||||
mu_of_r_dft_average = 0.d0
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
|
||||
rho = dm_a + dm_b
|
||||
if(mu_of_r_dft(i).gt.1.d+3)cycle
|
||||
mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i)
|
||||
enddo
|
||||
mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num)
|
||||
print*,'mu_of_r_dft_average = ',mu_of_r_dft_average
|
||||
END_PROVIDER
|
||||
|
||||
37
src/dft_utils_func/mu_of_r_dft.irp.f
Normal file
37
src/dft_utils_func/mu_of_r_dft.irp.f
Normal file
@ -0,0 +1,37 @@
|
||||
double precision function mu_rs_c(rho)
|
||||
implicit none
|
||||
double precision, intent(in) :: rho
|
||||
include 'constants.include.F'
|
||||
double precision :: cst_rs,alpha_rs,rs
|
||||
cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
|
||||
alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
|
||||
|
||||
rs = cst_rs * rho**(-1.d0/3.d0)
|
||||
mu_rs_c = alpha_rs/dsqrt(rs)
|
||||
|
||||
end
|
||||
|
||||
double precision function mu_grad_rho_func(r)
|
||||
implicit none
|
||||
double precision , intent(in) :: r(3)
|
||||
integer :: m
|
||||
double precision :: rho, dm_a, dm_b, grad_dm_a(3), grad_dm_b(3)
|
||||
double precision :: eta, grad_rho(3), grad_sqr
|
||||
eta = mu_erf
|
||||
call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
|
||||
rho = dm_a + dm_b
|
||||
do m = 1,3
|
||||
grad_rho(m) = grad_dm_a(m) + grad_dm_b(m)
|
||||
enddo
|
||||
grad_sqr=0.d0
|
||||
do m = 1,3
|
||||
grad_sqr=grad_sqr+grad_rho(m)*grad_rho(m)
|
||||
enddo
|
||||
grad_sqr = dsqrt(grad_sqr)
|
||||
if (rho<1.d-12) then
|
||||
mu_grad_rho_func = 1.d-10
|
||||
else
|
||||
mu_grad_rho_func = eta * grad_sqr / rho
|
||||
endif
|
||||
|
||||
end
|
||||
39
src/dft_utils_in_r/ints_grad.irp.f
Normal file
39
src/dft_utils_in_r/ints_grad.irp.f
Normal file
@ -0,0 +1,39 @@
|
||||
BEGIN_PROVIDER [ double precision, mo_grad_ints, (mo_num, mo_num,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mo_grad_ints(i,j,m) = <phi_i^MO | d/dx | phi_j^MO>
|
||||
END_DOC
|
||||
integer :: i,j,ipoint,m
|
||||
double precision :: weight
|
||||
mo_grad_ints = 0.d0
|
||||
do m = 1, 3
|
||||
do ipoint = 1, n_points_final_grid
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
do j = 1, mo_num
|
||||
do i = 1, mo_num
|
||||
mo_grad_ints(i,j,m) += mos_grad_in_r_array(j,ipoint,m) * mos_in_r_array(i,ipoint) * weight
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mo_grad_ints_transp, (3,mo_num, mo_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mo_grad_ints(i,j,m) = <phi_i^MO | d/dx | phi_j^MO>
|
||||
END_DOC
|
||||
integer :: i,j,ipoint,m
|
||||
double precision :: weight
|
||||
do m = 1, 3
|
||||
do j = 1, mo_num
|
||||
do i = 1, mo_num
|
||||
mo_grad_ints_transp(m,i,j) = mo_grad_ints(i,j,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
@ -21,7 +21,9 @@
|
||||
weight = final_weight_at_r_vector(i)
|
||||
rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate)
|
||||
rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate)
|
||||
call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
|
||||
energy_x_sr_lda(istate) += weight * e_x
|
||||
enddo
|
||||
enddo
|
||||
@ -48,7 +50,9 @@
|
||||
weight = final_weight_at_r_vector(i)
|
||||
rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate)
|
||||
rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate)
|
||||
call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
|
||||
energy_c_sr_lda(istate) += weight * e_c
|
||||
enddo
|
||||
enddo
|
||||
@ -122,8 +126,10 @@ END_PROVIDER
|
||||
weight = final_weight_at_r_vector(i)
|
||||
rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate)
|
||||
rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate)
|
||||
call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b)
|
||||
call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b)
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b)
|
||||
call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b)
|
||||
do j =1, ao_num
|
||||
aos_sr_vc_alpha_lda_w(j,i,istate) = sr_vc_a * aos_in_r_array(j,i)*weight
|
||||
aos_sr_vc_beta_lda_w(j,i,istate) = sr_vc_b * aos_in_r_array(j,i)*weight
|
||||
@ -147,8 +153,6 @@ END_PROVIDER
|
||||
double precision :: mu,weight
|
||||
double precision :: e_c,sr_vc_a,sr_vc_b,e_x,sr_vx_a,sr_vx_b
|
||||
double precision, allocatable :: rhoa(:),rhob(:)
|
||||
double precision :: mu_local
|
||||
mu_local = mu_erf_dft
|
||||
allocate(rhoa(N_states), rhob(N_states))
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
@ -158,6 +162,8 @@ END_PROVIDER
|
||||
weight = final_weight_at_r_vector(i)
|
||||
rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate)
|
||||
rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate)
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b)
|
||||
call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b)
|
||||
do j =1, ao_num
|
||||
|
||||
@ -36,8 +36,10 @@
|
||||
grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
|
||||
enddo
|
||||
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
! inputs
|
||||
call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
energy_x_sr_pbe(istate) += ex * weight
|
||||
@ -135,8 +137,10 @@ END_PROVIDER
|
||||
grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
|
||||
enddo
|
||||
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
! inputs
|
||||
call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
vx_rho_a *= weight
|
||||
@ -292,8 +296,10 @@ END_PROVIDER
|
||||
grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
|
||||
enddo
|
||||
|
||||
double precision :: mu_local
|
||||
mu_local = mu_of_r_dft(i)
|
||||
! inputs
|
||||
call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
vx_rho_a *= weight
|
||||
|
||||
@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
|
||||
enddo
|
||||
endif
|
||||
|
||||
call print_energy_components()
|
||||
! call print_energy_components()
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
@ -17,7 +17,7 @@ program rs_ks_scf
|
||||
print*, '**************************'
|
||||
print*, 'mu_erf_dft = ',mu_erf_dft
|
||||
print*, '**************************'
|
||||
call check_coherence_functional
|
||||
! call check_coherence_functional
|
||||
call create_guess
|
||||
call orthonormalize_mos
|
||||
call run
|
||||
|
||||
@ -1,6 +1,9 @@
|
||||
|
||||
subroutine give_all_mos_at_r(r,mos_array)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mos_array(i) = ith MO function evaluated at "r"
|
||||
END_DOC
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: mos_array(mo_num)
|
||||
double precision :: aos_array(ao_num)
|
||||
|
||||
@ -302,21 +302,21 @@ end
|
||||
integer(key_kind) :: idx
|
||||
double precision :: tmp
|
||||
|
||||
icount = 1 ! Avoid division by zero
|
||||
do j=1,mo_num
|
||||
do i=1,j-1
|
||||
call two_e_integrals_index(i,j,j,i,idx)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(mo_integrals_map,idx,tmp)
|
||||
banned_excitation(i,j) = dabs(tmp) < 1.d-14
|
||||
banned_excitation(j,i) = banned_excitation(i,j)
|
||||
if (banned_excitation(i,j)) icount = icount+2
|
||||
enddo
|
||||
enddo
|
||||
use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1%
|
||||
if (use_banned_excitation) then
|
||||
print *, 'Using sparsity of exchange integrals'
|
||||
endif
|
||||
!icount = 1 ! Avoid division by zero
|
||||
!do j=1,mo_num
|
||||
! do i=1,j-1
|
||||
! call two_e_integrals_index(i,j,j,i,idx)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call map_get(mo_integrals_map,idx,tmp)
|
||||
! banned_excitation(i,j) = dabs(tmp) < 1.d-14
|
||||
! banned_excitation(j,i) = banned_excitation(i,j)
|
||||
! if (banned_excitation(i,j)) icount = icount+2
|
||||
! enddo
|
||||
!enddo
|
||||
!use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1%
|
||||
!if (use_banned_excitation) then
|
||||
! print *, 'Using sparsity of exchange integrals'
|
||||
!endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -10,7 +10,7 @@ type:integer
|
||||
interface: ezfio,provider
|
||||
|
||||
[pseudo_n_k]
|
||||
doc: Number of gaussians in the local component
|
||||
doc: Powers of r - 2 in the local component
|
||||
type: integer
|
||||
interface: ezfio,provider
|
||||
size: (nuclei.nucl_num,pseudo.pseudo_klocmax)
|
||||
@ -38,7 +38,7 @@ type:integer
|
||||
interface: ezfio,provider
|
||||
|
||||
[pseudo_n_kl]
|
||||
doc: Number of functions in the non-local component
|
||||
doc: Powers of r - 2 in the non-local component
|
||||
type: integer
|
||||
interface: ezfio,provider
|
||||
size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax)
|
||||
|
||||
@ -2,3 +2,4 @@ fci
|
||||
mo_two_e_erf_ints
|
||||
aux_quantities
|
||||
hartree_fock
|
||||
two_body_rdm
|
||||
|
||||
@ -52,8 +52,8 @@ program molden
|
||||
l += 1
|
||||
if (l > ao_num) exit
|
||||
enddo
|
||||
write(i_unit_output,*)''
|
||||
enddo
|
||||
write(i_unit_output,*)''
|
||||
enddo
|
||||
|
||||
|
||||
|
||||
24
src/tools/save_natorb_no_ref.irp.f
Normal file
24
src/tools/save_natorb_no_ref.irp.f
Normal file
@ -0,0 +1,24 @@
|
||||
program save_natorb
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Save natural |MOs| into the |EZFIO|.
|
||||
!
|
||||
! This program reads the wave function stored in the |EZFIO| directory,
|
||||
! extracts the corresponding natural orbitals and setd them as the new
|
||||
! |MOs|.
|
||||
!
|
||||
! If this is a multi-state calculation, the density matrix that produces
|
||||
! the natural orbitals is obtained from an average of the density
|
||||
! matrices of each state with the corresponding
|
||||
! :option:`determinants state_average_weight`
|
||||
END_DOC
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
call save_natural_mos_canon_label
|
||||
call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None')
|
||||
call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None')
|
||||
call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None')
|
||||
call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None')
|
||||
call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None')
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user