From 31bb892b657fe7a2054cc5247fbe81fdc8e09978 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Jan 2024 11:45:24 +0100 Subject: [PATCH 01/12] Better error message --- ocaml/Zmatrix.ml | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/ocaml/Zmatrix.ml b/ocaml/Zmatrix.ml index 9e6ab2f8..6427f734 100644 --- a/ocaml/Zmatrix.ml +++ b/ocaml/Zmatrix.ml @@ -58,17 +58,32 @@ let int_of_atom_id : atom_id -> int = fun x -> x let float_of_distance : float StringMap.t -> distance -> float = fun map -> function | Value x -> x - | Label s -> StringMap.find s map + | Label s -> begin + try StringMap.find s map with + | Not_found -> + Printf.sprintf "Zmatrix error: distance %s undefined" s + |> failwith + end let float_of_angle : float StringMap.t -> angle -> float = fun map -> function | Value x -> x - | Label s -> StringMap.find s map + | Label s -> begin + try StringMap.find s map with + | Not_found -> + Printf.sprintf "Zmatrix error: angle %s undefined" s + |> failwith + end let float_of_dihedral : float StringMap.t -> dihedral -> float = fun map -> function | Value x -> x - | Label s -> StringMap.find s map + | Label s -> begin + try StringMap.find s map with + | Not_found -> + Printf.sprintf "Zmatrix error: dihedral %s undefined" s + |> failwith + end type line = From c0a4b7890e51454e078ab894b935a4e772484fab Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Jan 2024 13:19:21 +0100 Subject: [PATCH 02/12] Fix bug in complex svd --- src/utils/linear_algebra.irp.f | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 314ad4f6..7cef9ee4 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -645,7 +645,7 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff) END_DOC integer, intent(in) :: m,n, LDA, LDC complex*16, intent(in) :: A(LDA,n) - double precision, intent(in) :: cutoff + double precision, intent(in) :: cutoff, d1 complex*16, intent(out) :: C(LDC,m) double precision, allocatable :: D(:), rwork(:) @@ -673,8 +673,9 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff) stop 1 endif + d1 = D(1) do i=1,n - if (D(i) > cutoff*D(1)) then + if (D(i) > cutoff*d1) then D(i) = 1.d0/D(i) else D(i) = 0.d0 From 0b83c1ab8b34bd303142f0a7352b0775510ee874 Mon Sep 17 00:00:00 2001 From: ydamour Date: Fri, 26 Jan 2024 17:34:16 +0100 Subject: [PATCH 03/12] mkl with gfortran --- config/gfortran_mkl.cfg | 62 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 config/gfortran_mkl.cfg diff --git a/config/gfortran_mkl.cfg b/config/gfortran_mkl.cfg new file mode 100644 index 00000000..f2787d63 --- /dev/null +++ b/config/gfortran_mkl.cfg @@ -0,0 +1,62 @@ +# Common flags +############## +# +# -ffree-line-length-none : Needed for IRPF90 which produces long lines +# -lblas -llapack : Link with libblas and liblapack libraries provided by the system +# -I . : Include the curent directory (Mandatory) +# +# --ninja : Allow the utilisation of ninja. (Mandatory) +# --align=32 : Align all provided arrays on a 32-byte boundary +# +# +[COMMON] +FC : gfortran -ffree-line-length-none -I . -mavx -g -fPIC -std=legacy +LAPACK_LIB : -I${MKLROOT}/include -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_core -lpthread -lm -ldl -lmkl_gnu_thread -lgomp -fopenmp +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -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 +#################### +# +# -Ofast : Disregard strict standards compliance. Enables all -O3 optimizations. +# It also enables optimizations that are not valid +# for all standard-compliant programs. It turns on +# -ffast-math and the Fortran-specific +# -fno-protect-parens and -fstack-arrays. +[OPT] +FCFLAGS : -Ofast -mavx + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -Ofast + +# Debugging flags +################# +# +# -fcheck=all : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# +[DEBUG] +FCFLAGS : -fcheck=all -g + +# OpenMP flags +################# +# +[OPENMP] +FC : -fopenmp +IRPF90_FLAGS : --openmp + From b5b0cdb27a734162c2f0ab90e0aa83b33d13d490 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2024 08:50:14 +0100 Subject: [PATCH 04/12] README.md --- README.md | 4 ++++ external/ezfio | 2 +- external/irpf90 | 2 +- src/ao_one_e_ints/spread_dipole_ao.irp.f | 2 +- 4 files changed, 7 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index b03f2ecc..5a35f63d 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,7 @@ +**Important**: The Intel ifx compiler is not able to produce correct +executables for Quantum Package. Please use ifort as long as you can, and +consider switching to gfortran in the long term. + # Quantum Package 2.2 diff --git a/external/ezfio b/external/ezfio index d5805497..dba01c4f 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 +Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3 diff --git a/external/irpf90 b/external/irpf90 index 0007f72f..4ab1b175 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 +Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 diff --git a/src/ao_one_e_ints/spread_dipole_ao.irp.f b/src/ao_one_e_ints/spread_dipole_ao.irp.f index c52d0548..86469a3f 100644 --- a/src/ao_one_e_ints/spread_dipole_ao.irp.f +++ b/src/ao_one_e_ints/spread_dipole_ao.irp.f @@ -224,7 +224,7 @@ subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) BEGIN_DOC ! Computes the following integral : -! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x^2 ] ! needed for the dipole and those things END_DOC implicit none From 5b5df61960aad048452bd398c1ec584fffa9c267 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2024 14:13:10 +0100 Subject: [PATCH 05/12] Fixed linear algebra --- config/ifort_2021_avx.cfg | 2 +- config/ifort_2021_avx_mpi.cfg | 2 +- config/ifort_2021_avx_notz.cfg | 2 +- config/ifort_2021_debug.cfg | 2 +- config/ifort_2021_mpi_rome.cfg | 2 +- config/ifort_2021_rome.cfg | 2 +- config/ifort_2021_sse4.cfg | 2 +- config/ifort_2021_sse4_mpi.cfg | 2 +- config/ifort_2021_xHost.cfg | 2 +- src/utils/linear_algebra.irp.f | 3 ++- 10 files changed, 11 insertions(+), 10 deletions(-) diff --git a/config/ifort_2021_avx.cfg b/config/ifort_2021_avx.cfg index 6c34cf47..55fe0ee7 100644 --- a/config/ifort_2021_avx.cfg +++ b/config/ifort_2021_avx.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_avx_mpi.cfg b/config/ifort_2021_avx_mpi.cfg index 4c893c73..362f482a 100644 --- a/config/ifort_2021_avx_mpi.cfg +++ b/config/ifort_2021_avx_mpi.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : mpiifort -fpic +FC : mpiifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_avx_notz.cfg b/config/ifort_2021_avx_notz.cfg index 1fa595d7..3cd80236 100644 --- a/config/ifort_2021_avx_notz.cfg +++ b/config/ifort_2021_avx_notz.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 --define=WITHOUT_TRAILZ --define=WITHOUT_SHIFTRL diff --git a/config/ifort_2021_debug.cfg b/config/ifort_2021_debug.cfg index 80802f33..2e30642c 100644 --- a/config/ifort_2021_debug.cfg +++ b/config/ifort_2021_debug.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL diff --git a/config/ifort_2021_mpi_rome.cfg b/config/ifort_2021_mpi_rome.cfg index e47a466e..b7341388 100644 --- a/config/ifort_2021_mpi_rome.cfg +++ b/config/ifort_2021_mpi_rome.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : mpiifort -fpic +FC : mpiifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_rome.cfg b/config/ifort_2021_rome.cfg index 504438c9..1d2d8c77 100644 --- a/config/ifort_2021_rome.cfg +++ b/config/ifort_2021_rome.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_sse4.cfg b/config/ifort_2021_sse4.cfg index 07c3ebb8..e43147ba 100644 --- a/config/ifort_2021_sse4.cfg +++ b/config/ifort_2021_sse4.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic +FC : ifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_sse4_mpi.cfg b/config/ifort_2021_sse4_mpi.cfg index f3fa0eaa..1914988b 100644 --- a/config/ifort_2021_sse4_mpi.cfg +++ b/config/ifort_2021_sse4_mpi.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : mpiifort -fpic +FC : mpiifort -fpic -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_xHost.cfg b/config/ifort_2021_xHost.cfg index 9170b059..0dfce550 100644 --- a/config/ifort_2021_xHost.cfg +++ b/config/ifort_2021_xHost.cfg @@ -6,7 +6,7 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic -diag-disable 5462 +FC : ifort -fpic -diag-disable=5462 -diag-disable=10448 LAPACK_LIB : -mkl=parallel IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=64 -DINTEL diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 075525d1..c9d0be72 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -645,13 +645,14 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff) END_DOC integer, intent(in) :: m,n, LDA, LDC complex*16, intent(in) :: A(LDA,n) - double precision, intent(in) :: cutoff, d1 + double precision, intent(in) :: cutoff complex*16, intent(out) :: C(LDC,m) double precision, allocatable :: D(:), rwork(:) complex*16, allocatable :: U(:,:), Vt(:,:), work(:), A_tmp(:,:) integer :: info, lwork integer :: i,j,k + double precision :: d1 allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n),rwork(5*n)) do j=1,n do i=1,m From 419ed79c49a89382f3e473df647c9624a4f3e759 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 10 Feb 2024 12:48:29 +0100 Subject: [PATCH 06/12] added transition two rdm --- src/davidson/u0_wee_u0.irp.f | 22 + src/two_body_rdm/act_2_transition_rdm.irp.f | 39 + src/two_body_rdm/example.irp.f | 88 ++ src/two_body_rdm/io_two_rdm.irp.f | 34 + src/two_body_rdm/test_2_rdm.irp.f | 1 + .../davidson_like_trans_2rdm.irp.f | 585 ++++++++++ src/two_rdm_routines/update_trans_rdm.irp.f | 1002 +++++++++++++++++ 7 files changed, 1771 insertions(+) create mode 100644 src/two_body_rdm/act_2_transition_rdm.irp.f create mode 100644 src/two_rdm_routines/davidson_like_trans_2rdm.irp.f create mode 100644 src/two_rdm_routines/update_trans_rdm.irp.f diff --git a/src/davidson/u0_wee_u0.irp.f b/src/davidson/u0_wee_u0.irp.f index 0c543aca..bd3525e1 100644 --- a/src/davidson/u0_wee_u0.irp.f +++ b/src/davidson/u0_wee_u0.irp.f @@ -492,3 +492,25 @@ subroutine u_0_H_u_0_two_e(e_0,u_0,n,keys_tmp,Nint,N_st,sze) deallocate (s_0, v_0) end +BEGIN_PROVIDER [double precision, psi_energy_two_e_trans, (N_states, N_states)] + implicit none + BEGIN_DOC +! psi_energy_two_e_trans(istate,jstate) = + END_dOC + integer :: i,j,istate,jstate + double precision :: hij, coef_i, coef_j + psi_energy_two_e_trans = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_H_j_two_e(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + do istate = 1, N_states + coef_i = psi_coef(i,istate) + do jstate = 1, N_states + coef_j = psi_coef(j,jstate) + psi_energy_two_e_trans(jstate,istate) += coef_i * coef_j * hij + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/two_body_rdm/act_2_transition_rdm.irp.f b/src/two_body_rdm/act_2_transition_rdm.irp.f new file mode 100644 index 00000000..3d08b084 --- /dev/null +++ b/src/two_body_rdm/act_2_transition_rdm.irp.f @@ -0,0 +1,39 @@ + BEGIN_PROVIDER [double precision, act_2_rdm_trans_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states,N_states)] + implicit none + BEGIN_DOC +! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2rdm_trans +! +! \sum_{\sigma,\sigma'} +! +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1) +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for beta/beta spin + print*,'' + print*,'Providing act_2_rdm_trans_spin_trace_mo ' + character*(128) :: name_file + name_file = 'act_2_rdm_trans_spin_trace_mo' + ispin = 4 + act_2_rdm_trans_spin_trace_mo = 0.d0 + call wall_time(wall_1) +! if(read_two_body_rdm_trans_spin_trace)then +! print*,'Reading act_2_rdm_trans_spin_trace_mo from disk ...' +! call read_array_two_rdm_trans(n_act_orb,N_states,act_2_rdm_trans_spin_trace_mo,name_file) +! else + call orb_range_2_trans_rdm_openmp(act_2_rdm_trans_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) +! endif +! if(write_two_body_rdm_trans_spin_trace)then +! print*,'Writing act_2_rdm_trans_spin_trace_mo on disk ...' +! call write_array_two_rdm_trans(n_act_orb,n_states,act_2_rdm_trans_spin_trace_mo,name_file) +! call ezfio_set_two_body_rdm_trans_io_two_body_rdm_trans_spin_trace("Read") +! endif + + act_2_rdm_trans_spin_trace_mo *= 2.d0 + call wall_time(wall_2) + print*,'Wall time to provide act_2_rdm_trans_spin_trace_mo',wall_2 - wall_1 + END_PROVIDER diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 30e2685a..38510fe9 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -365,3 +365,91 @@ subroutine routine_full_mos end + +subroutine routine_active_only_trans + implicit none + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate,jstate + BEGIN_DOC +! This routine computes the two electron repulsion within the active space using various providers +! + END_DOC + + double precision :: vijkl,get_two_e_integral + double precision :: wee_tot(N_states,N_states),rdm_transtot + double precision :: spin_trace + double precision :: accu_tot + + wee_tot = 0.d0 + + + iorb = 1 + jorb = 1 + korb = 1 + lorb = 1 + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + provide act_2_rdm_trans_spin_trace_mo + i = 1 + j = 2 + + print*,'**************************' + print*,'**************************' + do jstate = 1, N_states + do istate = 1, N_states + !! PURE ACTIVE PART + !! + accu_tot = 0.d0 + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + ! 1 2 1 2 2 1 2 1 +! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate)).gt.1.d-10)then +! print*,'Error in act_2_rdm_trans_spin_trace_mo' +! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l) - act_2_rdm_trans_spin_trace_mo(j,i,l,k)).gt.1.d-10" +! print*,i,j,k,l +! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate)) +! endif + + ! 1 2 1 2 1 2 1 2 +! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10)then +! print*,'Error in act_2_rdm_trans_spin_trace_mo' +! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10" +! print*,i,j,k,l +! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)) +! endif + + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + + + rdm_transtot = act_2_rdm_trans_spin_trace_mo(l,k,j,i,istate,jstate) + + wee_tot(istate,jstate) += 0.5d0 * vijkl * rdm_transtot + + enddo + enddo + enddo + enddo + print*,'' + print*,'' + print*,'Active space only energy for state ',istate,jstate + print*,'wee_tot = ',wee_tot(istate,jstate) + print*,'Full energy ' + print*,'psi_energy_two_e(istate,jstate)= ',psi_energy_two_e_trans(istate,jstate) + print*,'--------------------------' + enddo + enddo + print*,'Wee from DM ' + do istate = 1,N_states + write(*,'(100(F16.10,X))')wee_tot(1:N_states,istate) + enddo + print*,'Wee from Psi det' + do istate = 1,N_states + write(*,'(100(F16.10,X))')psi_energy_two_e_trans(1:N_states,istate) + enddo + +end + diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f index bdd8a4f9..0b30d76f 100644 --- a/src/two_body_rdm/io_two_rdm.irp.f +++ b/src/two_body_rdm/io_two_rdm.irp.f @@ -31,3 +31,37 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file) close(unit=i_unit_output) end + +subroutine write_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates,nstates) + + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'W') + call lock_io() + write(i_unit_output)array_tmp + call unlock_io() + close(unit=i_unit_output) +end + +subroutine read_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states,nstates) + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'R') + call lock_io() + read(i_unit_output)array_tmp + call unlock_io() + close(unit=i_unit_output) +end + diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index 123261d8..de2606a7 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -4,5 +4,6 @@ program test_2_rdm touch read_wf call routine_active_only call routine_full_mos + call routine_active_only_trans end diff --git a/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f new file mode 100644 index 00000000..9e68a0e1 --- /dev/null +++ b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f @@ -0,0 +1,585 @@ +subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sze) + use bitmasks + implicit none + BEGIN_DOC + ! if ispin == 1 :: alpha/alpha 2_rdm + ! == 2 :: beta /beta 2_rdm + ! == 3 :: alpha/beta + beta/alpha 2trans_rdm + ! == 4 :: spin traced 2_rdm :: aa + bb + ab + ba + ! + ! notice that here it is the TRANSITION RDM THAT IS COMPUTED + ! + ! THE DIAGONAL PART IS THE USUAL ONE FOR A GIVEN STATE + ! Assumes that the determinants are in psi_det + ! + ! istart, iend, ishift, istep are used in ZMQ parallelization. + END_DOC + integer, intent(in) :: N_st,sze + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) + double precision, intent(in) :: u_0(sze,N_st) + + integer :: k + double precision, allocatable :: u_t(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t + PROVIDE mo_two_e_integrals_in_map + allocate(u_t(N_st,N_det)) + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) + enddo + call dtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_st) + + call orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1) + deallocate(u_t) + + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + enddo + +end + +subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes two-trans_rdm + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) + double precision, intent(in) :: u_t(N_st,N_det) + + integer :: k + + PROVIDE N_int + + select case (N_int) + case (1) + call orb_range_2_trans_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (2) + call orb_range_2_trans_rdm_openmp_work_2(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (3) + call orb_range_2_trans_rdm_openmp_work_3(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (4) + call orb_range_2_trans_rdm_openmp_work_4(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case default + call orb_range_2_trans_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + end select +end + + + BEGIN_TEMPLATE +subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + use omp_lib + implicit none + BEGIN_DOC + ! Computes the two trans_rdm for the N_st vectors |u_t> + ! if ispin == 1 :: alpha/alpha 2trans_rdm + ! == 2 :: beta /beta 2trans_rdm + ! == 3 :: alpha/beta 2trans_rdm + ! == 4 :: spin traced 2trans_rdm :: aa + bb + 0.5 (ab + ba)) + ! The 2trans_rdm will be computed only on the list of orbitals list_orb, which contains norb + ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + double precision, intent(in) :: u_t(N_st,N_det) + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) + + integer(omp_lock_kind) :: lock_2trans_rdm + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b + integer :: krow, kcol + integer :: lrow, lcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + integer(bit_kind) :: orb_bitmask($N_int) + integer :: list_orb_reverse(mo_num) + integer, allocatable :: keys(:,:) + double precision, allocatable :: values(:,:,:) + integer :: nkeys,sze_buff + integer :: ll + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + else + print*,'Wrong parameter for ispin in general_2_trans_rdm_state_av_openmp_work' + print*,'ispin = ',ispin + stop + endif + + + PROVIDE N_int + + call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) + sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + list_orb_reverse = -1000 + do i = 1, norb + list_orb_reverse(list_orb(i)) = i + enddo + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + call omp_init_lock(lock_2trans_rdm) + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + PROVIDE N_int nthreads_davidson elec_alpha_num + !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & + !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2trans_rdm,& + !$OMP psi_bilinear_matrix_columns, & + !$OMP psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,& + !$OMP psi_bilinear_matrix_transp_rows, & + !$OMP psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_transp_order, N_st, & + !$OMP psi_bilinear_matrix_order_transp_reverse, & + !$OMP psi_bilinear_matrix_columns_loc, & + !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & + !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & + !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & + !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, & + !$OMP lcol, lrow, l_a, l_b, & + !$OMP buffer, doubles, n_doubles, & + !$OMP tmp_det2, idx, l, kcol_prev, & + !$OMP singles_a, n_singles_a, singles_b, & + !$OMP n_singles_b, nkeys, keys, values) + + ! Alpha/Beta double excitations + ! ============================= + nkeys = 0 + allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff)) + allocate( buffer($N_int,maxab), & + singles_a(maxab), & + singles_b(maxab), & + doubles(maxab), & + idx(maxab)) + + kcol_prev=-1 + + ASSERT (iend <= N_det) + ASSERT (istart > 0) + ASSERT (istep > 0) + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + if (kcol /= kcol_prev) then + call get_all_spin_singles_$N_int( & + psi_det_beta_unique, idx0, & + tmp_det(1,2), N_det_beta_unique, & + singles_b, n_singles_b) + endif + kcol_prev = kcol + + ! Loop over singly excited beta columns + ! ------------------------------------- + + do i=1,n_singles_b + lcol = singles_b(i) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) + + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) + + ASSERT (l_a <= N_det) + idx(j) = l_a + l_a = l_a+1 + enddo + j = j-1 + + call get_all_spin_singles_$N_int( & + buffer, idx, tmp_det(1,1), j, & + singles_a, n_singles_a ) + + ! Loop over alpha singles + ! ----------------------- + + if(alpha_beta.or.spin_trace)then + do k = 1,n_singles_a + l_a = singles_a(k) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) +! print*,'nkeys before = ',nkeys + do ll = 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta)then + ! only ONE contribution + if (nkeys+1 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + else if (spin_trace)then + ! TWO contributions + if (nkeys+2 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + endif + call orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + + enddo + endif + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + enddo + + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha exitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta.or.spin_trace.or.alpha_alpha)then + ! increment the alpha/beta part for single excitations + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the alpha/alpha part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + + enddo + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + if(alpha_alpha.or.spin_trace)then + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + enddo + endif + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + l_a = psi_bilinear_matrix_transp_order(l_b) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta.or.spin_trace.or.beta_beta)then + ! increment the alpha/beta part for single excitations + if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the beta /beta part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_bb_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + enddo + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + if(beta_beta.or.spin_trace)then + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + l_a = psi_bilinear_matrix_transp_order(l_b) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +! print*,'to do orb_range_off_diag_double_to_2_trans_rdm_bb_dm_buffer' + ASSERT (l_a <= N_det) + + enddo + endif + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem + + double precision :: c_1(N_states,N_states) + do ll = 1, N_states + do l = 1, N_states + c_1(l,ll) = u_t(ll,k_a) * u_t(l,k_a) + enddo + enddo + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + call orb_range_diag_to_all_states_2_rdm_trans_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + end do + !$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) + !$OMP END PARALLEL + +end + + SUBST [ N_int ] + + 1;; + 2;; + 3;; + 4;; + N_int;; + + END_TEMPLATE + +subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + use omp_lib + implicit none + integer, intent(in) :: n_st,nkeys,dim1 + integer, intent(in) :: keys(4,nkeys) + double precision, intent(in) :: values(n_st,n_st,nkeys) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st,n_st) + + integer(omp_lock_kind),intent(inout):: lock_2rdm + + integer :: i,h1,h2,p1,p2,istate,jstate + call omp_set_lock(lock_2rdm) + +! print*,'*************' +! print*,'updating' +! print*,'nkeys',nkeys + do i = 1, nkeys + h1 = keys(1,i) + h2 = keys(2,i) + p1 = keys(3,i) + p2 = keys(4,i) + do jstate = 1, N_st + do istate = 1, N_st +!! print*,h1,h2,p1,p2,values(istate,i) + big_array(h1,h2,p1,p2,istate,jstate) += values(istate,jstate,i) + enddo + enddo + enddo + call omp_unset_lock(lock_2rdm) + +end + diff --git a/src/two_rdm_routines/update_trans_rdm.irp.f b/src/two_rdm_routines/update_trans_rdm.irp.f new file mode 100644 index 00000000..9f7077a2 --- /dev/null +++ b/src/two_rdm_routines/update_trans_rdm.irp.f @@ -0,0 +1,1002 @@ + subroutine orb_range_diag_to_all_states_2_rdm_trans_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for a given determinant det_1 + ! + ! c_1 is the array of the contributions to the trans_rdm for all states + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: det_1(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2 + integer(bit_kind) :: det_1_act(N_int,2) + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + do i = 1, N_int + det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + enddo + + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) + logical :: is_integer_in_string + integer :: i1,i2,istate + integer :: jstate + if(alpha_beta)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + ! If alpha/beta, electron 1 is alpha, electron 2 is beta + ! Therefore you don't necessayr have symmetry between electron 1 and 2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + + else if (alpha_alpha)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if (beta_beta)then + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if(spin_trace)then + ! 0.5 * (alpha beta + beta alpha) + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + endif + end + + + subroutine orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another +! +! c_1 is the array of the contributions to the trans_rdm for all states +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-trans_rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-trans_rdm +! +! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + integer :: jstate + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation(det_1,det_2,exc,phase,N_int) + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 = exc(1,1,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 = exc(1,2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_beta)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + else if(spin_trace)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + endif + end + + subroutine orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body trans_rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if(alpha_beta)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + endif + else if(spin_trace)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + endif + endif + end + + subroutine orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body trans_rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + use bitmasks + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if(alpha_alpha.or.spin_trace)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + else + return + endif + endif + end + + subroutine orb_range_off_diag_single_to_all_states_bb_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body trans_rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if(beta_beta.or.spin_trace)then + if (exc(0,1,1) == 1) then + return + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + endif + endif + end + + + subroutine orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body trans_rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + + integer :: jstate + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_alpha.or.spin_trace)then + nkeys += 1 + + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end + + subroutine orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body trans_rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(beta_beta.or.spin_trace)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end + From 1b9a75f4886a4112920da3c4f611e19bf35cae12 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 12 Feb 2024 18:18:53 +0100 Subject: [PATCH 07/12] Fixed pseudo-inverse (extrapolations) --- src/mol_properties/EZFIO.cfg | 7 +++ .../print_e_components.irp.f | 0 src/mol_properties/print_mol_properties.irp.f | 7 ++- src/utils/linear_algebra.irp.f | 44 +++++++++---------- 4 files changed, 34 insertions(+), 24 deletions(-) rename src/{two_body_rdm => mol_properties}/print_e_components.irp.f (100%) diff --git a/src/mol_properties/EZFIO.cfg b/src/mol_properties/EZFIO.cfg index 35a095fb..3ddba227 100644 --- a/src/mol_properties/EZFIO.cfg +++ b/src/mol_properties/EZFIO.cfg @@ -21,3 +21,10 @@ type: logical doc: If true and N_states > 1, the oscillator strength will be computed interface: ezfio,provider,ocaml default: false + +[calc_energy_components] +type: logical +doc: If true, the components of the energy (1e, 2e, kinetic) will be computed +interface: ezfio,provider,ocaml +default: false + diff --git a/src/two_body_rdm/print_e_components.irp.f b/src/mol_properties/print_e_components.irp.f similarity index 100% rename from src/two_body_rdm/print_e_components.irp.f rename to src/mol_properties/print_e_components.irp.f diff --git a/src/mol_properties/print_mol_properties.irp.f b/src/mol_properties/print_mol_properties.irp.f index 3753a3dd..00ccb826 100644 --- a/src/mol_properties/print_mol_properties.irp.f +++ b/src/mol_properties/print_mol_properties.irp.f @@ -6,6 +6,11 @@ subroutine print_mol_properties() ! Run the propertie calculations END_DOC + ! Energy components + if (calc_energy_components) then + call print_energy_components + endif + ! Electric dipole moment if (calc_dipole_moment) then call print_dipole_moment @@ -18,7 +23,7 @@ subroutine print_mol_properties() ! Oscillator strength if (calc_osc_str .and. N_states > 1) then - call print_oscillator_strength + call print_oscillator_strength endif end diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index c9d0be72..76b280b7 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1377,31 +1377,29 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff) enddo endif - print*, ' n_svd = ', n_svd - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j) & - !$OMP SHARED (n, n_svd, D, Vt) - !$OMP DO - do j = 1, n - do i = 1, n_svd - Vt(i,j) = D(i) * Vt(i,j) - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - call dgemm("N", "N", m, n, n_svd, 1.d0, U, m, Vt, n, 0.d0, C, LDC) - -! C = 0.d0 -! do i=1,m -! do j=1,n -! do k=1,n -! C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) -! enddo +! !$OMP PARALLEL & +! !$OMP DEFAULT (NONE) & +! !$OMP PRIVATE (i, j) & +! !$OMP SHARED (n, n_svd, D, Vt) +! !$OMP DO +! do j = 1, n +! do i = 1, n_svd +! Vt(i,j) = D(i) * Vt(i,j) ! enddo ! enddo +! !$OMP END DO +! !$OMP END PARALLEL + +! call dgemm('N', 'N', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1)) + + C = 0.d0 + do i=1,m + do j=1,n + do k=1,n_svd + C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) + enddo + enddo + enddo deallocate(U,D,Vt,work,A_tmp) From d619c621fcd00e35d8bb8c2f956486044366a6d0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 12 Feb 2024 18:21:59 +0100 Subject: [PATCH 08/12] DGEMM in pseudo-inverse --- src/utils/linear_algebra.irp.f | 42 +++++++++++++++++----------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 76b280b7..2db47092 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1377,29 +1377,29 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff) enddo endif -! !$OMP PARALLEL & -! !$OMP DEFAULT (NONE) & -! !$OMP PRIVATE (i, j) & -! !$OMP SHARED (n, n_svd, D, Vt) -! !$OMP DO -! do j = 1, n -! do i = 1, n_svd -! Vt(i,j) = D(i) * Vt(i,j) -! enddo -! enddo -! !$OMP END DO -! !$OMP END PARALLEL - -! call dgemm('N', 'N', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1)) - - C = 0.d0 - do i=1,m - do j=1,n - do k=1,n_svd - C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) - enddo + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j) & + !$OMP SHARED (n, n_svd, D, Vt) + !$OMP DO + do j = 1, n + do i = 1, n_svd + Vt(i,j) = D(i) * Vt(i,j) enddo enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm('T', 'T', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1)) + +! C = 0.d0 +! do i=1,m +! do j=1,n +! do k=1,n_svd +! C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) +! enddo +! enddo +! enddo deallocate(U,D,Vt,work,A_tmp) From fbb946d8f44161bb8de5c752060e82980265717b Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 15 Feb 2024 16:46:05 +0100 Subject: [PATCH 09/12] removed the systematic save of MOs in casscf --- src/casscf_cipsi/casscf.irp.f | 2 +- src/two_body_rdm/act_2_transition_rdm.irp.f | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index addca236..c0cd361d 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -99,8 +99,8 @@ subroutine run mo_coef = NewOrbs mo_occ = occnum - call save_mos if(.not.converged)then + call save_mos iteration += 1 if(norm_grad_vec2.gt.0.01d0)then N_det = N_states diff --git a/src/two_body_rdm/act_2_transition_rdm.irp.f b/src/two_body_rdm/act_2_transition_rdm.irp.f index 3d08b084..612213e2 100644 --- a/src/two_body_rdm/act_2_transition_rdm.irp.f +++ b/src/two_body_rdm/act_2_transition_rdm.irp.f @@ -1,9 +1,9 @@ BEGIN_PROVIDER [double precision, act_2_rdm_trans_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states,N_states)] implicit none BEGIN_DOC -! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2rdm_trans +! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) = STATE SPECIFIC physicist notation for 2rdm_trans ! -! \sum_{\sigma,\sigma'} +! \sum_{\sigma,\sigma'} ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! From 22c99a0484eb75ed85c789fa7e39bc965c7fd591 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 15 Feb 2024 19:32:15 +0100 Subject: [PATCH 10/12] done some cleaning in the casscf and added a detailed example of Multi state CASSCF --- src/casscf_cipsi/README.rst | 9 ++- src/casscf_cipsi/casscf.irp.f | 16 ++++- src/casscf_cipsi/example_casscf_multistate.sh | 66 +++++++++++++++++++ 3 files changed, 85 insertions(+), 6 deletions(-) create mode 100755 src/casscf_cipsi/example_casscf_multistate.sh diff --git a/src/casscf_cipsi/README.rst b/src/casscf_cipsi/README.rst index f84cde75..75c99de2 100644 --- a/src/casscf_cipsi/README.rst +++ b/src/casscf_cipsi/README.rst @@ -4,13 +4,15 @@ casscf |CASSCF| program with the CIPSI algorithm. -Example of inputs ------------------ + +Example of inputs for GROUND STATE calculations +----------------------------------------------- +NOTICE :: FOR EXCITED STATES CALCULATIONS SEE THE FILE "example_casscf_multistate.sh" a) Small active space : standard CASSCF --------------------------------------- Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units) -3 +2 O 0.0000000000 0.0000000000 -1.1408000000 O 0.0000000000 0.0000000000 1.1408000000 @@ -45,3 +47,4 @@ qp set casscf_cipsi small_active_space False qp run casscf | tee ${EZFIO_FILE}.casscf_large.out # you should find around -149.9046 + diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index c0cd361d..d0a26d36 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -54,14 +54,24 @@ subroutine run call write_time(6) call write_int(6,iteration,'CAS-SCF iteration = ') - call write_double(6,energy,'CAS-SCF energy = ') + call write_double(6,energy,'State-average CAS-SCF energy = ') ! if(n_states == 1)then ! call ezfio_get_casscf_cipsi_energy_pt2(E_PT2) ! call ezfio_get_casscf_cipsi_energy(PT2) + double precision :: delta_E_istate, e_av + e_av = 0.d0 do istate=1,N_states - call write_double(6,E_PT2(istate),'E + PT2 energy = ') - call write_double(6,PT2(istate),' PT2 = ') + e_av += state_average_weight(istate) * Ev(istate) + if(istate.gt.1)then + delta_E_istate = E_PT2(istate) - E_PT2(1) + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2 = ',delta_E_istate + endif + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate) + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' PT2 energy = ',PT2(istate) +! call write_double(6,E_PT2(istate),'E + PT2 energy = ') +! call write_double(6,PT2(istate),' PT2 = ') enddo + call write_double(6,e_av,'State-average CAS-SCF energy bis = ') call write_double(6,pt2_max,' PT2_MAX = ') ! endif diff --git a/src/casscf_cipsi/example_casscf_multistate.sh b/src/casscf_cipsi/example_casscf_multistate.sh new file mode 100755 index 00000000..368c0440 --- /dev/null +++ b/src/casscf_cipsi/example_casscf_multistate.sh @@ -0,0 +1,66 @@ +# This is an example for MULTI STATE CALCULATION STATE AVERAGE CASSCF +# We will compute 3 states on the O2 molecule +# The Ground state and 2 degenerate excited states +# Please follow carefully the tuto :) + +##### PREPARING THE EZFIO +# Set the path to your QP2 directory +QP_ROOT=my_fancy_path +source ${QP_ROOT}/quantum_package.rc +# Create the EZFIO folder +qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz_multi_state +# Start with ROHF orbitals +qp run scf +# Freeze the 1s orbitals of the two oxygen +qp set_frozen_core + +##### PREPARING THE ORBITALS WITH NATURAL ORBITALS OF A CIS +# Tell that you want 3 states in your WF +qp set determinants n_states 3 +# Run a CIS wave function to start your calculation +qp run cis | tee ${EZFIO_FILE}.cis_3_states.out +# Save the STATE AVERAGE natural orbitals for having a balanced description +# This will also order the orbitals according to their occupation number +# Which makes the active space selection easyer ! +qp run save_natorb | tee ${EZFIO_FILE}.natorb_3states.out + +##### PREPARING A CIS GUESS WITHIN THE ACTIVE SPACE +# Set an active space which has the most of important excitations +# and that maintains symmetry : the ACTIVE ORBITALS are from """6 to 13""" + +# YOU FIRST FREEZE THE VIRTUALS THAT ARE NOT IN THE ACTIVE SPACE +# !!!!! WE SET TO "-D" for DELETED !!!! +qp set_mo_class -c "[1-5]" -a "[6-13]" -d "[14-46]" +# You create a guess of CIS type WITHIN THE ACTIVE SPACE +qp run cis | tee ${EZFIO_FILE}.cis_3_states_active_space.out +# You tell to read the WFT stored (i.e. the guess we just created) +qp set determinants read_wf True + +##### DOING THE CASSCF +### SETTING PROPERLY THE ACTIVE SPACE FOR CASSCF +# You set the active space WITH THE VIRTUAL ORBITALS !!! +# !!!!! NOW WE SET TO "-v" for VIRTUALS !!!!! +qp set_mo_class -c "[1-5]" -a "[6-13]" -v "[14-46]" + +# You tell that it is a small actice space so the CIPSI can take all Slater determinants +qp set casscf_cipsi small_active_space True +# You specify the output file +output=${EZFIO_FILE}.casscf_3states.out +# You run the CASSCF calculation +qp run casscf | tee ${output} + +# Some grep in order to get some numbers useful to check convergence +# State average energy +grep "State-average CAS-SCF energy =" $output | cut -d "=" -f 2 > data_e_average +# Delta E anticipated for State-average energy, only usefull to check convergence +grep "Predicted energy improvement =" $output | cut -d "=" -f 2 > data_improve +# Ground state energy +grep "state 1 E + PT2 energy" $output | cut -d "=" -f 2 > data_1 +# First excited state energy +grep "state 2 E + PT2 energy" $output | cut -d "=" -f 2 > data_2 +# First excitation energy +grep "state 2 Delta E+PT2" $output | cut -d "=" -f 2 > data_delta_E2 +# Second excited state energy +grep "state 3 E + PT2 energy" $output | cut -d "=" -f 2 > data_3 +# Second excitation energy +grep "state 3 Delta E+PT2" $output | cut -d "=" -f 2 > data_delta_E3 From fa877df399a918750c28a0a262d27823c0cbd3c6 Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Feb 2024 15:12:39 +0100 Subject: [PATCH 11/12] added exponential of anti-hermitian matrices using the Helgaker's book formulation, and of general matrices using the Taylor expansion. Replaced in casscf_cipsi Umat variable --- src/casscf_cipsi/neworbs.irp.f | 41 +++++----- src/utils/linear_algebra.irp.f | 137 +++++++++++++++++++++++++++++++++ 2 files changed, 158 insertions(+), 20 deletions(-) diff --git a/src/casscf_cipsi/neworbs.irp.f b/src/casscf_cipsi/neworbs.irp.f index a7cebbb2..ca2deebb 100644 --- a/src/casscf_cipsi/neworbs.irp.f +++ b/src/casscf_cipsi/neworbs.irp.f @@ -226,27 +226,28 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] end do ! Form the exponential + call exp_matrix_taylor(Tmat,mo_num,Umat,converged) - Tpotmat(:,:)=0.D0 - Umat(:,:) =0.D0 - do i=1,mo_num - Tpotmat(i,i)=1.D0 - Umat(i,i) =1.d0 - end do - iter=0 - converged=.false. - do while (.not.converged) - iter+=1 - f = 1.d0 / dble(iter) - Tpotmat2(:,:) = Tpotmat(:,:) * f - call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, & - Tpotmat2, size(Tpotmat2,1), & - Tmat, size(Tmat,1), 0.d0, & - Tpotmat, size(Tpotmat,1)) - Umat(:,:) = Umat(:,:) + Tpotmat(:,:) - - converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) - end do +! Tpotmat(:,:)=0.D0 +! Umat(:,:) =0.D0 +! do i=1,mo_num +! Tpotmat(i,i)=1.D0 +! Umat(i,i) =1.d0 +! end do +! iter=0 +! converged=.false. +! do while (.not.converged) +! iter+=1 +! f = 1.d0 / dble(iter) +! Tpotmat2(:,:) = Tpotmat(:,:) * f +! call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, & +! Tpotmat2, size(Tpotmat2,1), & +! Tmat, size(Tmat,1), 0.d0, & +! Tpotmat, size(Tpotmat,1)) +! Umat(:,:) = Umat(:,:) + Tpotmat(:,:) +! +! converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) +! end do END_PROVIDER diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 2db47092..175beff3 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1897,3 +1897,140 @@ end do end subroutine pivoted_cholesky +subroutine exp_matrix(X,n,exp_X) + implicit none + double precision, intent(in) :: X(n,n) + integer, intent(in):: n + double precision, intent(out):: exp_X(n,n) + BEGIN_DOC + ! exponential of the matrix X: X has to be ANTI HERMITIAN !! + ! + ! taken from Hellgaker, jorgensen, Olsen book + ! + ! section evaluation of matrix exponential (Eqs. 3.1.29 to 3.1.31) + END_DOC + integer :: i + double precision, allocatable :: r2_mat(:,:),eigvalues(:),eigvectors(:,:) + double precision, allocatable :: matrix_tmp1(:,:),eigvalues_mat(:,:),matrix_tmp2(:,:) + include 'constants.include.F' + allocate(r2_mat(n,n),eigvalues(n),eigvectors(n,n)) + allocate(eigvalues_mat(n,n),matrix_tmp1(n,n),matrix_tmp2(n,n)) + + ! r2_mat = X^2 in the 3.1.30 + call get_A_squared(X,n,r2_mat) + call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n) + eigvalues=-eigvalues + + if(.False.)then + !!! For debugging and following the book intermediate + ! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30 + ! matrix_tmp1 = W t^2 + print*,'eigvalues = ' + do i = 1, n + print*,i,eigvalues(i) + write(*,'(100(F16.10,X))')eigvectors(:,i) + enddo + eigvalues_mat=0.d0 + do i = 1,n + ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 + eigvalues(i) = dsqrt(eigvalues(i)) + eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i) + enddo + call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & + eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) + call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & + eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) + print*,'r2_mat new = ' + do i = 1, n + write(*,'(100(F16.10,X))')matrix_tmp2(:,i) + enddo + endif + + ! building the exponential + ! exp(X) = W cos(t) W^T + W t^-1 sin(t) W^T X as in Eq. 3.1.31 + ! matrix_tmp1 = W cos(t) + do i = 1,n + eigvalues_mat(i,i) = dcos(eigvalues(i)) + enddo + ! matrix_tmp2 = W cos(t) + call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & + eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) + ! matrix_tmp2 = W cos(t) W^T + call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & + eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) + exp_X = matrix_tmp2 + ! matrix_tmp2 = W t^-1 sin(t) W^T X + do i = 1,n + if(dabs(eigvalues(i)).gt.1.d-4)then + eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i) + else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6 + eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 + endif + enddo + ! matrix_tmp1 = W t^-1 sin(t) + call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & + eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) + ! matrix_tmp2 = W t^-1 sin(t) W^T + call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & + eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) + ! exp_X += matrix_tmp2 X + call dgemm('N','N',n,n,n,1.d0,matrix_tmp2,size(matrix_tmp2,1), & + X,size(X,1),1.d0,exp_X,size(exp_X,1)) + +end + + +subroutine exp_matrix_taylor(X,n,exp_X,converged) + implicit none + BEGIN_DOC + ! exponential of a general real matrix X using the Taylor expansion of exp(X) + ! + ! returns the logical converged which checks the convergence + END_DOC + double precision, intent(in) :: X(n,n) + integer, intent(in):: n + double precision, intent(out):: exp_X(n,n) + logical :: converged + double precision :: f + integer :: i,iter + double precision, allocatable :: Tpotmat(:,:),Tpotmat2(:,:) + allocate(Tpotmat(n,n),Tpotmat2(n,n)) + BEGIN_DOC + ! exponential of X using Taylor expansion + END_DOC + Tpotmat(:,:)=0.D0 + exp_X(:,:) =0.D0 + do i=1,n + Tpotmat(i,i)=1.D0 + exp_X(i,i) =1.d0 + end do + iter=0 + converged=.false. + do while (.not.converged) + iter+=1 + f = 1.d0 / dble(iter) + Tpotmat2(:,:) = Tpotmat(:,:) * f + call dgemm('N','N', n,n,n,1.d0, & + Tpotmat2, size(Tpotmat2,1), & + X, size(X,1), 0.d0, & + Tpotmat, size(Tpotmat,1)) + exp_X(:,:) = exp_X(:,:) + Tpotmat(:,:) + + converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) + end do + if(.not.converged)then + print*,'Warning !! exp_matrix_taylor did not converge !' + endif + +end + +subroutine get_A_squared(A,n,A2) + implicit none + BEGIN_DOC +! A2 = A A where A is n x n matrix. Use the dgemm routine + END_DOC + double precision, intent(in) :: A(n,n) + integer, intent(in) :: n + double precision, intent(out):: A2(n,n) + call dgemm('N','N',n,n,n,1.d0,A,size(A,1),A,size(A,1),0.d0,A2,size(A2,1)) +end From ac805f9f016ab5b035557642e19602144d845c6f Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Feb 2024 15:25:38 +0100 Subject: [PATCH 12/12] added some reference numbers in the example_casscf_multistate.sh --- src/casscf_cipsi/example_casscf_multistate.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/casscf_cipsi/example_casscf_multistate.sh b/src/casscf_cipsi/example_casscf_multistate.sh index 368c0440..716c211a 100755 --- a/src/casscf_cipsi/example_casscf_multistate.sh +++ b/src/casscf_cipsi/example_casscf_multistate.sh @@ -10,7 +10,7 @@ source ${QP_ROOT}/quantum_package.rc # Create the EZFIO folder qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz_multi_state # Start with ROHF orbitals -qp run scf +qp run scf # ROHF energy : -149.619992871398 # Freeze the 1s orbitals of the two oxygen qp set_frozen_core @@ -18,7 +18,7 @@ qp set_frozen_core # Tell that you want 3 states in your WF qp set determinants n_states 3 # Run a CIS wave function to start your calculation -qp run cis | tee ${EZFIO_FILE}.cis_3_states.out +qp run cis | tee ${EZFIO_FILE}.cis_3_states.out # -149.6652601409258 -149.4714726176746 -149.4686165431939 # Save the STATE AVERAGE natural orbitals for having a balanced description # This will also order the orbitals according to their occupation number # Which makes the active space selection easyer ! @@ -32,7 +32,7 @@ qp run save_natorb | tee ${EZFIO_FILE}.natorb_3states.out # !!!!! WE SET TO "-D" for DELETED !!!! qp set_mo_class -c "[1-5]" -a "[6-13]" -d "[14-46]" # You create a guess of CIS type WITHIN THE ACTIVE SPACE -qp run cis | tee ${EZFIO_FILE}.cis_3_states_active_space.out +qp run cis | tee ${EZFIO_FILE}.cis_3_states_active_space.out # -149.6515472533511 -149.4622878024821 -149.4622878024817 # You tell to read the WFT stored (i.e. the guess we just created) qp set determinants read_wf True @@ -47,7 +47,7 @@ qp set casscf_cipsi small_active_space True # You specify the output file output=${EZFIO_FILE}.casscf_3states.out # You run the CASSCF calculation -qp run casscf | tee ${output} +qp run casscf | tee ${output} # -149.7175867510 -149.5059010227 -149.5059010226 # Some grep in order to get some numbers useful to check convergence # State average energy