From 57eabff6758b254d9c6e92b04e356205509ecc1d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 29 Jun 2019 17:29:32 +0200 Subject: [PATCH] added documentation for the two-rdm --- src/casscf/densities.irp.f | 8 +- src/casscf/get_energy.irp.f | 14 +- src/two_body_rdm/ab_only_routines.irp.f | 22 +-- src/two_body_rdm/all_2rdm_routines.irp.f | 2 +- src/two_body_rdm/orb_range_2_rdm.irp.f | 20 +++ ...outines.irp.f => orb_range_routines.irp.f} | 0 src/two_body_rdm/routines_compute_2rdm.irp.f | 14 +- .../routines_compute_2rdm_orb_range.irp.f | 129 +++++++++++++++++- 8 files changed, 168 insertions(+), 41 deletions(-) rename src/two_body_rdm/{general_2rdm_routines.irp.f => orb_range_routines.irp.f} (100%) diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 7b243bb4..30a914f1 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -42,7 +42,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 if (bavard) then - write(6,*) ' providing density matrix P0' + write(6,*) ' providing the 2 body RDM on the active part' endif P0tuvx= 0.d0 @@ -55,11 +55,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] uu = list_act(u) do t = 1, n_act_orb tt = list_act(t) - P0tuvx(t,u,v,x) = & - state_average_weight(istate) * & - ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) + & - two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) + & - two_rdm_beta_beta_mo (tt,uu,vv,xx,istate) ) + P0tuvx(t,u,v,x) = act_two_rdm_spin_trace_mo(t,v,u,x) enddo enddo enddo diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 29a12cad..0a5cfb49 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -1,5 +1,10 @@ program print_2rdm implicit none + BEGIN_DOC + ! get the active part of the bielectronic energy on a given wave function. + ! + ! useful to test the active part of the spin trace 2 rdms + END_DOC read_wf = .True. touch read_wf call routine @@ -23,18 +28,9 @@ subroutine routine i = list_act(ii) integral = get_two_e_integral(i,j,k,l,mo_integrals_map) accu(1) += act_two_rdm_spin_trace_mo(ii,jj,kk,ll) * integral - !if(dabs(act_two_rdm_spin_trace_mo(ii,jj,kk,ll)).gt.thr)then - !print*,'',ii,jj,kk,ll,act_two_rdm_spin_trace_mo(ii,jj,kk,ll)*integral - !print*,'accu',accu(1) - !endif enddo enddo enddo enddo print*,'accu = ',accu(1) - print*,'psi_energy_two_e = ',psi_energy_two_e -!double precision :: hij -!call i_H_j_double_alpha_beta(psi_det(1,1,1),psi_det(1,1,2),N_int,hij) -!print*,'hij * 2',hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 -!print*,'psi diag = ',psi_energy_two_e - hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 end diff --git a/src/two_body_rdm/ab_only_routines.irp.f b/src/two_body_rdm/ab_only_routines.irp.f index 195f439a..9041c753 100644 --- a/src/two_body_rdm/ab_only_routines.irp.f +++ b/src/two_body_rdm/ab_only_routines.irp.f @@ -1,9 +1,9 @@ - subroutine two_rdm_dm_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) + subroutine two_rdm_ab_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! Computes the alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS ! ! Assumes that the determinants are in psi_det ! @@ -27,7 +27,7 @@ size(u_t, 1), & N_det, N_st) - call two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) + call two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -37,11 +37,11 @@ end - subroutine two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! Computes the alpha/beta part of the two-body density matrix ! ! Default should be 1,N_det,0,1 END_DOC @@ -55,20 +55,20 @@ select case (N_int) case (1) - call two_rdm_dm_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call two_rdm_dm_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call two_rdm_dm_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call two_rdm_dm_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case default - call two_rdm_dm_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) end select end BEGIN_TEMPLATE - subroutine two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none integer, intent(in) :: N_st,sze,istart,iend,ishift,istep diff --git a/src/two_body_rdm/all_2rdm_routines.irp.f b/src/two_body_rdm/all_2rdm_routines.irp.f index 75d71ded..3f08b18f 100644 --- a/src/two_body_rdm/all_2rdm_routines.irp.f +++ b/src/two_body_rdm/all_2rdm_routines.irp.f @@ -2,7 +2,7 @@ subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab, use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! Computes the alpha/alpha, beta/beta and alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS ! ! Assumes that the determinants are in psi_det ! diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index e98612c5..c40c46d2 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -4,6 +4,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -17,6 +21,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -30,6 +38,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -48,6 +60,14 @@ BEGIN_PROVIDER [double precision, act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none + BEGIN_DOC +! act_two_rdm_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices +! The active part of the two-electron energy can be computed as: +! +! \sum_{i,j,k,l = 1, n_act_orb} act_two_rdm_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > +! +! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) + END_DOC double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) diff --git a/src/two_body_rdm/general_2rdm_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f similarity index 100% rename from src/two_body_rdm/general_2rdm_routines.irp.f rename to src/two_body_rdm/orb_range_routines.irp.f diff --git a/src/two_body_rdm/routines_compute_2rdm.irp.f b/src/two_body_rdm/routines_compute_2rdm.irp.f index 7165576f..112d2e36 100644 --- a/src/two_body_rdm/routines_compute_2rdm.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm.irp.f @@ -3,7 +3,7 @@ subroutine diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the DIAGONAL PART of the alpha/beta two body rdm +! routine that update the DIAGONAL PART of the alpha/beta two body rdm IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 @@ -31,7 +31,7 @@ subroutine diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the DIAGONAL PART of ALL THREE two body rdm +! routine that update the DIAGONAL PART of ALL THREE two body rdm IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 @@ -77,7 +77,7 @@ subroutine off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 @@ -101,7 +101,7 @@ subroutine off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 @@ -140,7 +140,7 @@ subroutine off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS END_DOC use bitmasks implicit none @@ -177,7 +177,7 @@ subroutine off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 @@ -214,7 +214,7 @@ subroutine off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS END_DOC implicit none integer, intent(in) :: dim1,dim2,dim3,dim4 diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f index c2283fb2..a3c7a76d 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -28,7 +28,20 @@ subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the DIAGONAL PART of ALL THREE two body rdm +! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1 +! +! 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -154,7 +167,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 3 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -219,7 +249,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 3 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -320,7 +367,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 1 or 4 will do something END_DOC use bitmasks implicit none @@ -383,7 +447,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 2 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -449,7 +530,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 1 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -505,7 +603,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body 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 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-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 2 or 4 will do something END_DOC implicit none