LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-13 16:15:18 +02:00
Code Releases Activity

Merge branch 'dev-stable' into dev-stable-tc-scf

Browse Source
This commit is contained in:
Abdallah Ammar 2023-03-04 14:33:43 +01:00
commit 4a23b63a42
34 changed files with 761 additions and 332 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
INSTALL.rst
View File

@ -316,7 +316,7 @@ OCaml
.. code:: bash
opam install ocamlbuild cryptokit zmq sexplib ppx_sexp_conv ppx_deriving getopt
opam install ocamlbuild zmq sexplib ppx_sexp_conv ppx_deriving getopt
Docopt

66
codemeta.json Normal file
View File

@ -0,0 +1,66 @@
{
"@context": "https://doi.org/10.5063/schema/codemeta-2.0",
"@type": "SoftwareSourceCode",
"license": "https://spdx.org/licenses/AGPL-3.0",
"codeRepository": "https://github.com/QuantumPackage/qp2",
"dateCreated": "2014-04-01",
"datePublished": "2019-06-11",
"dateModified": "2020-02-21",
"downloadUrl": "https://github.com/QuantumPackage/qp2/releases/tag/2.1.2",
"issueTracker": "https://github.com/QuantumPackage/qp2/issues",
"name": "Quantum Package",
"version": "2.1.2",
"identifier": "https://doi.org/10.5281/zenodo.3677565",
"description": "Programming environment for wave function methods",
"applicationCategory": "Quantum Chemistry",
"funding": "ERC_863481, CoE_952165",
"developmentStatus": "active",
"referencePublication": "https://doi.org/10.1021/acs.jctc.9b00176",
"funder": {
"@type": "Organization",
"name": "CNRS"
},
"keywords": [
"selected configuration interaction",
"CIPSI"
],
"programmingLanguage": [
"Fortran",
"IRPF90",
"OCaml",
"Python",
"C"
],
"operatingSystem": [
"Linux"
],
"softwareRequirements": [
"ZeroMQ"
],
"author": [
{
"@type": "Person",
"@id": "https://orcid.org/0000-0003-4955-7136",
"givenName": "Scemama",
"familyName": "Anthony",
"email": "scemama@irsamc.ups-tlse.fr",
"affiliation": {
"@type": "Organization",
"name": "Laboratoire de chimie et physique quantiques, Toulouse, CNRS"
}
}
],
"contributor": [
{
"@type": "Person",
"givenName": "Emmanuel",
"familyName": "Giner",
"email": "eginer@lct.jussieu.fr",
"affiliation": {
"@type": "Organization",
"name": "Laboratoire de Chimie Theorique, Paris, CNRS"
}
}
]
}

2
configure vendored
View File

@ -19,6 +19,8 @@ git submodule update
# Update ARM or x86 dependencies
ARCHITECTURE=$(uname -m)
cd ${QP_ROOT}/external/qp2-dependencies
git checkout master
git pull
echo "Architecture: $ARCHITECTURE"
cd ${QP_ROOT}

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit b8cd5815bce14c9b880e3c5ea3d5fc2652f5955c
Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8

2
ocaml/Makefile
View File

@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m
$(QP_ROOT)/ocaml/element_create_db.byte
external_libs:
opam install cryptokit sexplib
opam install sexplib
qpackage.odocl: $(MLIFILES)
ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl

4
ocaml/To_md5.ml
View File

@ -4,8 +4,8 @@ open Sexplib
let to_md5 sexp_of_t t =
sexp_of_t t
|> Sexp.to_string
|> Cryptokit.hash_string (Cryptokit.Hash.md5 ())
|> Cryptokit.transform_string (Cryptokit.Hexa.encode ())
|> Digest.string
|> Digest.to_hex
|> MD5.of_string
;;

2
ocaml/_tags
View File

@ -1,4 +1,4 @@
true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
true: package(zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
true: thread
false: profile
<*byte> : linkdep(c_bindings.o), custom

3
src/basis_correction/TODO
View File

@ -1 +1,4 @@
change all correlation functionals with the pbe_on_top general
factor 2 in two-rdm involved in:
on-top, mu(r), pbe-on-top, sc_basis_corr and so on
casscf : state_av_act_2_rdm_spin_trace_mo

14
src/basis_correction/pbe_on_top.irp.f
View File

@ -41,13 +41,13 @@
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
! We take the extrapolated on-top pair density
on_top_extrap = mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
@ -108,8 +108,8 @@
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
! We take the extrapolated on-top pair density
on_top_extrap = mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
@ -170,8 +170,8 @@
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We DO NOT take the extrapolated on-top pair density, but there is * 2 because of normalization
on_top_extrap = 2.d0 * on_top
! We DO NOT take the extrapolated on-top pair density
on_top_extrap = on_top
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)

6
src/basis_correction/print_routine.irp.f
View File

@ -18,7 +18,7 @@ subroutine print_basis_correction
print*, ''
print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018) '
print*, ' Journal of Physical Chemistry Letters 10, 2931-2937 (2019) '
print*, ' ???REF SC?'
print*, ' Journal of Chemical Physics 152, 174104 (2020) '
print*, '****************************************'
print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential
@ -56,14 +56,14 @@ subroutine print_basis_correction
print*,''
print*,'********************************************'
print*,'********************************************'
print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)'
print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)'
print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)

6
src/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -46,9 +46,9 @@ END_PROVIDER
&BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_y , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_z , (mo_num,mo_num)]
BEGIN_DOC
! array of the integrals of MO_i * x MO_j
! array of the integrals of MO_i * y MO_j
! array of the integrals of MO_i * z MO_j
! array of the integrals of Left MO_i * x Right MO_j
! array of the integrals of Left MO_i * y Right MO_j
! array of the integrals of Left MO_i * z Right MO_j
END_DOC
implicit none

2
src/dft_utils_in_r/dm_in_r.irp.f
View File

@ -82,7 +82,7 @@
enddo
enddo
!$OMP END PARALLEL DO
print*,'density and gradients provided'
! print*,'density and gradients provided'
END_PROVIDER

20
src/hartree_fock/10.hf.bats
View File

@ -34,23 +34,28 @@ cat > hcn_charges.xyz << EOF
0.5 -2.0 0.0 0.0
EOF
rm -rf hcn.ezfio
qp create_ezfio -b def2-svp hcn.xyz
EZFIO=hcn_pt_charges
rm -rf $EZFIO
qp create_ezfio -b def2-svp hcn.xyz -o $EZFIO
qp run scf
mv hcn_charges.xyz hcn.ezfio_point_charges.xyz
python write_pt_charges.py hcn.ezfio
mv hcn_charges.xyz ${EZFIO}_point_charges.xyz
python write_pt_charges.py ${EZFIO}
qp set nuclei point_charges True
qp run scf | tee hcn.ezfio.pt_charges.out
qp run scf | tee ${EZFIO}.pt_charges.out
energy="$(ezfio get hartree_fock energy)"
rm -rf hcn.ezfio
good=-92.76613324421798
eq $energy $good $thresh
rm -rf $EZFIO
}
@test "point charges" {
run_pt_charges
}
@test "HCN" { # 7.792500 8.51926s
run hcn.ezfio -92.88717500035233
}
@test "B-B" { # 3s
run b2_stretched.ezfio -48.9950585434279
}
@ -124,9 +129,6 @@ good=-92.76613324421798
run ch4.ezfio -40.19961807784367
}
@test "HCN" { # 7.792500 8.51926s
run hcn.ezfio -92.88717500035233
}
@test "N2" { # 8.648100 13.754s
run n2.ezfio -108.9834897852979

5
src/mu_of_r/f_hf_utils.irp.f
View File

@ -86,6 +86,9 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
enddo
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
f_HF_val_ab *= 2.d0
two_bod_dens *= 2.d0
end
@ -136,4 +139,6 @@ subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab)
enddo
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
int_f_HF_val_ab *= 2.d0
end

8
src/mu_of_r/f_psi_i_a_v_utils.irp.f
View File

@ -49,6 +49,9 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
enddo
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
f_ii_val_ab *= 2.d0
two_bod_dens *= 2.d0
end
@ -142,6 +145,9 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
f_ia_val_ab += v_tilde(i,a) * rho_tilde(i,a)
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
f_ia_val_ab *= 2.d0
two_bod_dens *= 2.d0
end
@ -222,6 +228,8 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
f_aa_val_ab += v_tilde(b,a) * rho_tilde(b,a)
enddo
enddo
! DO NOT multiply by two as in give_f_ii_val_ab and give_f_ia_val_ab because the N(N-1) normalization condition of the active two-rdm
end
BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)]

2
src/mu_of_r/f_psi_utils.irp.f
View File

@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_
do j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb
! 1 2 1 2
full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += 0.5d0 * full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo

2
src/mu_of_r/mu_of_r_conditions.irp.f
View File

@ -110,7 +110,7 @@
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
f_psi = f_psi_cas_ab(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
w_psi = 1.d+10
else

80
src/nuclei/write_pt_charges.py Normal file
View File

@ -0,0 +1,80 @@
#!/usr/bin/env python
import os
import sys
# First argument is the EZFIO file
# It reads a file EZFIO_point_charges.xyz written in this way:
# charge x y z (Angstrom)
# for all charges
def zip_in_ezfio(ezfio,tmp):
tmpzip=tmp+".gz"
cmdzip="gzip -c "+tmp+" > "+tmpzip
os.system(cmdzip)
os.system("rm "+tmp)
cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip
os.system(cmdmv)
def mv_in_ezfio(ezfio,tmp):
cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp
os.system(cmdmv)
# Getting the EZFIO
EZFIO=sys.argv[1]
EZFIO=EZFIO.replace("/", "")
print(EZFIO)
# Reading the point charges and convert the Angstrom geometry in Bohr for QP
f = open(EZFIO+'_point_charges.xyz','r')
lines = f.readlines()
convert_angs_to_bohr=1.8897259885789233
n_charges=0
coord_x=[]
coord_y=[]
coord_z=[]
charges=[]
for line in lines:
data = line.split()
if(len(data)>0):
n_charges += 1
charges.append(str(data[0]))
coord_x.append(str(convert_angs_to_bohr*float(data[1])))
coord_y.append(str(convert_angs_to_bohr*float(data[2])))
coord_z.append(str(convert_angs_to_bohr*float(data[3])))
# Write the file containing the number of charges and set in EZFIO folder
tmp="n_pts_charge"
fncharges = open(tmp,'w')
fncharges.write(" "+str(n_charges)+'\n')
fncharges.close()
mv_in_ezfio(EZFIO,tmp)
# Write the file containing the charges and set in EZFIO folder
tmp="pts_charge_z"
fcharges = open(tmp,'w')
fcharges.write(" 1\n")
fcharges.write(" "+str(n_charges)+'\n')
for i in range(n_charges):
fcharges.write(charges[i]+'\n')
fcharges.close()
zip_in_ezfio(EZFIO,tmp)
# Write the file containing the charge coordinates and set in EZFIO folder
tmp="pts_charge_coord"
fcoord = open(tmp,'w')
fcoord.write(" 2\n")
fcoord.write(" "+str(n_charges)+' 3\n')
#fcoord.write(" "+' 3 '+str(n_charges)+' \n')
for i in range(n_charges):
fcoord.write(' '+coord_x[i]+'\n')
for i in range(n_charges):
fcoord.write(' '+coord_y[i]+'\n')
for i in range(n_charges):
fcoord.write(' '+coord_z[i]+'\n')
fcoord.close()
zip_in_ezfio(EZFIO,tmp)

0
src/tc_bi_ortho/12.tc_bi_ortho.bats → src/tc_bi_ortho/31.tc_bi_ortho.bats
View File

16
src/tc_bi_ortho/print_tc_spin_dens.irp.f Normal file
View File

@ -0,0 +1,16 @@
program test_spin_dens
implicit none
BEGIN_DOC
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end.
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call tc_print_mulliken_sd
! call test
end

114
src/tc_bi_ortho/spin_mulliken.irp.f Normal file
View File

@ -0,0 +1,114 @@
BEGIN_PROVIDER [double precision, tc_spin_population, (ao_num,ao_num,N_states)]
implicit none
integer :: i,j,istate
BEGIN_DOC
! spin population on the ao basis :
! tc_spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
tc_spin_population = 0.d0
if(only_spin_tc_right)then
do i = 1, ao_num
do j = 1, ao_num
tc_spin_population(j,i,1) = tc_spin_dens_right_only(j,i) * ao_overlap(j,i)
enddo
enddo
else
do istate = 1, N_states
do i = 1, ao_num
do j = 1, ao_num
tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i)
enddo
enddo
enddo
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum, (0:ao_l_max,N_states)]
&BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum_per_atom, (0:ao_l_max,nucl_num,N_states)]
implicit none
integer :: i,istate
double precision :: accu
tc_spin_population_angular_momentum = 0.d0
tc_spin_population_angular_momentum_per_atom = 0.d0
do istate = 1, N_states
do i = 1, ao_num
tc_spin_population_angular_momentum(ao_l(i),istate) += tc_spin_gross_orbital_product(i,istate)
tc_spin_population_angular_momentum_per_atom(ao_l(i),ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_spin_gross_orbital_product, (ao_num,N_states)]
implicit none
tc_spin_gross_orbital_product = 0.d0
integer :: i,j,istate
BEGIN_DOC
! gross orbital product for the spin population
END_DOC
do istate = 1, N_states
do i = 1, ao_num
do j = 1, ao_num
tc_spin_gross_orbital_product(i,istate) += tc_spin_population(j,i,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_mulliken_spin_densities, (nucl_num,N_states)]
implicit none
integer :: i,j,istate
BEGIN_DOC
!ATOMIC SPIN POPULATION (ALPHA MINUS BETA)
END_DOC
tc_mulliken_spin_densities = 0.d0
do istate = 1, N_states
do i = 1, ao_num
tc_mulliken_spin_densities(ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate)
enddo
enddo
END_PROVIDER
subroutine tc_print_mulliken_sd
implicit none
double precision :: accu
integer :: i
integer :: j
print*,'Mulliken spin densities'
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),tc_mulliken_spin_densities(i,1)
accu += tc_mulliken_spin_densities(i,1)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += tc_spin_gross_orbital_product(i,1)
write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_character(ao_l(i))),tc_spin_gross_orbital_product(i,1)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += tc_spin_population_angular_momentum(i,1)
print*,' ',trim(l_to_character(i)),tc_spin_population_angular_momentum(i,1)
print*,'sum = ',accu
enddo
print*,'Angular momentum analysis per atom'
print*,'Angular momentum analysis'
do j = 1,nucl_num
accu = 0.d0
do i = 0, ao_l_max
accu += tc_spin_population_angular_momentum_per_atom(i,j,1)
write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_character(i)),tc_spin_population_angular_momentum_per_atom(i,j,1)
print*,'sum = ',accu
enddo
enddo
end

2
src/tc_bi_ortho/tc_natorb.irp.f
View File

@ -21,7 +21,7 @@
allocate(dm_tmp(mo_num,mo_num), fock_diag(mo_num))
dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1)
dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1)
print *, ' dm_tmp'
do i = 1, mo_num

145
src/tc_bi_ortho/tc_prop.irp.f
View File

@ -1,8 +1,11 @@
BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ]
BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_beta, (mo_num, mo_num,N_states,N_states) ]
&BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_alpha, (mo_num, mo_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_transition_matrix(p,h,istate,jstate) = <Chi_istate| a^\dagger_p a_h |Phi_jstate>
! tc_transition_matrix_mo_alpha(p,h,istate,jstate) = <Chi_istate| a^\dagger_p,alpha a_h,alpha |Phi_jstate>
!
! tc_transition_matrix_mo_beta(p,h,istate,jstate) = <Chi_istate| a^\dagger_p,beta a_h,beta |Phi_jstate>
!
! where <Chi_istate| and |Phi_jstate> are the left/right eigenvectors on a bi-ortho basis
END_DOC
@ -11,43 +14,65 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state
integer, allocatable :: occ(:,:)
integer :: n_occ_ab(2),degree,exc(0:2,2,2)
allocate(occ(N_int*bit_kind_size,2))
tc_transition_matrix = 0.d0
do istate = 1, N_states
do jstate = 1, N_states
tc_transition_matrix_mo_alpha = 0.d0
tc_transition_matrix_mo_beta = 0.d0
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree.gt.1)then
cycle
else if (degree == 0)then
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do p = 1, n_occ_ab(1) ! browsing the alpha electrons
m = occ(p,1)
tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
do p = 1, n_occ_ab(2) ! browsing the beta electrons
m = occ(p,1)
tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
else
call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Single alpha
h = exc(1,1,1) ! hole in psi_det(1,1,j)
p = exc(1,2,1) ! particle in psi_det(1,1,j)
else
! Single beta
h = exc(1,1,2) ! hole in psi_det(1,1,j)
p = exc(1,2,2) ! particle in psi_det(1,1,j)
endif
tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
endif
if(degree.gt.1)cycle
do istate = 1, N_states
do jstate = 1, N_states
if (degree == 0)then
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do p = 1, n_occ_ab(1) ! browsing the alpha electrons
m = occ(p,1)
tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
do p = 1, n_occ_ab(2) ! browsing the beta electrons
m = occ(p,1)
tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
else
call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Single alpha
h = exc(1,1,1) ! hole in psi_det(1,1,j)
p = exc(1,2,1) ! particle in psi_det(1,1,j)
tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
else
! Single beta
h = exc(1,1,2) ! hole in psi_det(1,1,j)
p = exc(1,2,2) ! particle in psi_det(1,1,j)
tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
endif
endif
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_transition_matrix_mo(p,h,istate,jstate) = \sum_{sigma=alpha,beta} <Chi_istate| a^\dagger_p,sigma a_h,sigma |Phi_jstate>
!
! where <Chi_istate| and |Phi_jstate> are the left/right eigenvectors on a bi-ortho basis
END_DOC
tc_transition_matrix_mo = tc_transition_matrix_mo_beta + tc_transition_matrix_mo_alpha
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_spin_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta
!
! where <Chi_istate| and |Phi_jstate> are the left/right eigenvectors on a bi-ortho basis
END_DOC
tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)]
implicit none
@ -57,9 +82,9 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i)
tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i)
tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i)
tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i)
tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i)
tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i)
enddo
enddo
enddo
@ -78,3 +103,55 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, tc_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_transition_matrix(p,h,istate,jstate) in the AO basis
END_DOC
integer :: i,j,k,l
double precision :: dm_mo
tc_transition_matrix_ao = 0.d0
integer :: istate,jstate
do istate = 1, N_states
do jstate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
dm_mo = tc_transition_matrix_mo(j,i,jstate,istate)
do k = 1, ao_num
do l = 1, ao_num
tc_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, tc_spin_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_spin_transition_matrix_ao(p,h,istate,jstate) in the AO basis
END_DOC
integer :: i,j,k,l
double precision :: dm_mo
tc_spin_transition_matrix_ao = 0.d0
integer :: istate,jstate
do istate = 1, N_states
do jstate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
dm_mo = tc_spin_transition_matrix_mo(j,i,jstate,istate)
do k = 1, ao_num
do l = 1, ao_num
tc_spin_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER

30
src/tc_bi_ortho/test_spin_dens.irp.f Normal file
View File

@ -0,0 +1,30 @@
BEGIN_PROVIDER [ double precision, mo_r_coef_normalized, (ao_num,mo_num) ]
implicit none
integer :: i,j,k
double precision :: norm
do i = 1, mo_num
norm = 0.d0
do j = 1, ao_num
do k = 1, ao_num
norm += mo_r_coef(k,i) * mo_r_coef(j,i) * ao_overlap(k,j)
enddo
enddo
norm = 1.d0/dsqrt(norm)
do j = 1, ao_num
mo_r_coef_normalized(j,i) = mo_r_coef(j,i) * norm
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, tc_spin_dens_right_only, (ao_num, ao_num)]
implicit none
integer :: i,j,k
tc_spin_dens_right_only = 0.d0
do i = elec_beta_num+1, elec_alpha_num
do j = 1, ao_num
do k = 1, ao_num
tc_spin_dens_right_only(k,j) += mo_r_coef_normalized(k,i) * mo_r_coef_normalized(j,i)
enddo
enddo
enddo
END_PROVIDER

6
src/tc_keywords/EZFIO.cfg
View File

@ -201,3 +201,9 @@ type: integer
doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
interface: ezfio,provider,ocaml
default: -1
[only_spin_tc_right]
type: logical
doc: If |true|, only the right part of WF is used to compute spin dens
interface: ezfio,provider,ocaml
default: False

23
src/two_body_rdm/act_2_rdm.irp.f
View File

@ -4,21 +4,18 @@
BEGIN_DOC
! 12 12
! 1 2 1 2 == <ij|kl>
! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons
! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons
!
! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
!
! + <Psi_{istate}| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi_{istate}>
!
! 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_{\alpha}^{act} * N_{\beta}^{act}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
!
! act_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
END_DOC
integer :: ispin
double precision :: wall_1, wall_2
@ -44,6 +41,7 @@
endif
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1
act_2_rdm_ab_mo *= 2.d0
END_PROVIDER
@ -56,7 +54,7 @@
!
! 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_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -84,6 +82,7 @@
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1
act_2_rdm_aa_mo *= 2.d0
END_PROVIDER
@ -96,7 +95,7 @@
!
! 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_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -124,6 +123,7 @@
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1
act_2_rdm_bb_mo *= 2.d0
END_PROVIDER
BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
@ -135,7 +135,7 @@
!
! 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)/2
! 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
@ -161,6 +161,7 @@
call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
endif
act_2_rdm_spin_trace_mo *= 2.d0
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1
END_PROVIDER

52
src/two_body_rdm/example.irp.f
View File

@ -117,10 +117,10 @@ subroutine routine_active_only
endif
wee_ab(istate) += vijkl * rdmab
wee_aa(istate) += vijkl * rdmaa
wee_bb(istate) += vijkl * rdmbb
wee_tot(istate) += vijkl * rdmtot
wee_ab(istate) += 0.5d0 * vijkl * rdmab
wee_aa(istate) += 0.5d0 * vijkl * rdmaa
wee_bb(istate) += 0.5d0 * vijkl * rdmbb
wee_tot(istate) += 0.5d0 * vijkl * rdmtot
enddo
enddo
@ -144,13 +144,13 @@ subroutine routine_active_only
print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate)
print*,'--------------------------'
print*,'accu_aa = ',accu_aa
print*,'N_a (N_a-1)/2 = ', elec_alpha_num*(elec_alpha_num-1)*0.5
print*,'N_a (N_a-1) = ', elec_alpha_num*(elec_alpha_num-1)
print*,'accu_bb = ',accu_bb
print*,'N_b (N_b-1)/2 = ', elec_beta_num*(elec_beta_num-1)*0.5
print*,'2 N_b (N_b-1) = ', elec_beta_num*(elec_beta_num-1)*2
print*,'accu_ab = ',accu_ab
print*,'N_a N_b = ', elec_beta_num*elec_alpha_num
print*,'accu_tot = ',accu_tot
print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num * 0.5
print*,'Ne(Ne-1) = ',(elec_num-1)*elec_num
enddo
wee_aa_st_av = 0.d0
wee_bb_st_av = 0.d0
@ -192,10 +192,10 @@ subroutine routine_active_only
print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av)
endif
wee_aa_st_av += vijkl * rdm_aa_st_av
wee_bb_st_av += vijkl * rdm_bb_st_av
wee_ab_st_av += vijkl * rdm_ab_st_av
wee_tot_st_av += vijkl * rdm_tot_st_av
wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av
wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av
wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av
wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
enddo
enddo
enddo
@ -279,10 +279,10 @@ subroutine routine_full_mos
rdmbb = full_occ_2_rdm_bb_mo(l,k,j,i,istate)
rdmtot = full_occ_2_rdm_spin_trace_mo(l,k,j,i,istate)
wee_ab(istate) += vijkl * rdmab
wee_aa(istate) += vijkl * rdmaa
wee_bb(istate) += vijkl * rdmbb
wee_tot(istate)+= vijkl * rdmtot
wee_ab(istate) += 0.5d0 * vijkl * rdmab
wee_aa(istate) += 0.5d0 * vijkl * rdmaa
wee_bb(istate) += 0.5d0 * vijkl * rdmbb
wee_tot(istate)+= 0.5d0 * vijkl * rdmtot
enddo
enddo
aa_norm(istate) += full_occ_2_rdm_aa_mo(j,i,j,i,istate)
@ -310,18 +310,19 @@ subroutine routine_full_mos
print*,'Normalization of two-rdms '
print*,''
print*,'aa_norm(istate) = ',aa_norm(istate)
print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(1) * (elec_num_tab(1) - 1)/2
print*,'N_alpha(N_alpha-1) = ',elec_num_tab(1) * (elec_num_tab(1) - 1)
print*,''
print*,'bb_norm(istate) = ',bb_norm(istate)
print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(2) * (elec_num_tab(2) - 1)/2
print*,'N_alpha(N_alpha-1) = ',elec_num_tab(2) * (elec_num_tab(2) - 1)
print*,''
print*,'ab_norm(istate) = ',ab_norm(istate)
print*,'N_alpha * N_beta = ',elec_num_tab(1) * elec_num_tab(2)
print*,'N_alpha * N_beta *2 = ',elec_num_tab(1) * elec_num_tab(2) * 2
print*,''
print*,'tot_norm(istate) = ',tot_norm(istate)
print*,'N(N-1)/2 = ',elec_num*(elec_num - 1)/2
print*,'N(N-1) = ',elec_num*(elec_num - 1)
enddo
! return
wee_aa_st_av = 0.d0
wee_bb_st_av = 0.d0
wee_ab_st_av = 0.d0
@ -341,10 +342,10 @@ subroutine routine_full_mos
rdm_ab_st_av = state_av_full_occ_2_rdm_ab_mo(l,k,j,i)
rdm_tot_st_av = state_av_full_occ_2_rdm_spin_trace_mo(l,k,j,i)
wee_aa_st_av += vijkl * rdm_aa_st_av
wee_bb_st_av += vijkl * rdm_bb_st_av
wee_ab_st_av += vijkl * rdm_ab_st_av
wee_tot_st_av += vijkl * rdm_tot_st_av
wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av
wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av
wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av
wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
enddo
enddo
enddo
@ -354,17 +355,12 @@ subroutine routine_full_mos
print*,''
print*,'STATE AVERAGE ENERGY '
print*,'wee_aa_st_av = ',wee_aa_st_av
print*,'wee_aa_st_av_2 = ',wee_aa_st_av_2
print*,'wee_bb_st_av = ',wee_bb_st_av
print*,'wee_bb_st_av_2 = ',wee_bb_st_av_2
print*,'wee_ab_st_av = ',wee_ab_st_av
print*,'wee_ab_st_av_2 = ',wee_ab_st_av_2
print*,'Sum of components = ',wee_aa_st_av + wee_bb_st_av + wee_ab_st_av
print*,'Sum of components_2 = ',wee_aa_st_av_2 + wee_bb_st_av_2 + wee_ab_st_av_2
print*,''
print*,'Full energy '
print*,'wee_tot_st_av = ',wee_tot_st_av
print*,'wee_tot_st_av_2 = ',wee_tot_st_av_2
print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3
end

154
src/two_body_rdm/full_orb_2_rdm.irp.f
View File

@ -4,22 +4,18 @@
full_occ_2_rdm_ab_mo = 0.d0
integer :: i,j,k,l,iorb,jorb,korb,lorb,istate
BEGIN_DOC
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta + beta/alpha electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! + <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA
!
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO ARE SET TO ZERO
END_DOC
full_occ_2_rdm_ab_mo = 0.d0
@ -50,7 +46,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -64,7 +60,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -76,7 +72,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
enddo
enddo
@ -93,7 +89,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -107,7 +103,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -119,7 +115,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
enddo
enddo
endif
@ -139,7 +135,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -172,11 +168,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -186,8 +182,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
@ -203,11 +199,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -217,8 +213,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -237,7 +233,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -270,11 +266,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -284,8 +280,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
@ -301,11 +297,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -315,8 +311,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -335,14 +331,14 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
! The two-electron energy of each state can be computed as:
!
! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll >
! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * 1/2 * < ii jj | kk ll >
!
! with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l)
END_DOC
@ -377,11 +373,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -390,8 +386,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -403,11 +399,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -416,8 +412,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -433,11 +429,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -446,8 +442,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -459,11 +455,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -472,8 +468,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -489,14 +485,14 @@
korb = list_inact(k)
! ALPHA INACTIVE - BETA ACTIVE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! BETA INACTIVE - ALPHA ACTIVE
! beta alph beta alpha
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -506,8 +502,8 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
enddo
enddo
@ -523,14 +519,14 @@
korb = list_core(k)
!! BETA ACTIVE - ALPHA CORE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
!! ALPHA ACTIVE - BETA CORE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -540,8 +536,8 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
enddo
enddo

21
src/two_body_rdm/state_av_act_2rdm.irp.f
View File

@ -2,21 +2,16 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha/beta+beta/alpha electron pairs
!
! = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
!
! 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_{\alpha}^{act} * N_{\beta}^{act}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
!
! state_av_act_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
@ -34,6 +29,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_ab_mo *= 2.d0
END_PROVIDER
@ -48,7 +45,7 @@
!
! 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_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -63,6 +60,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_aa_mo *= 2.d0
END_PROVIDER
@ -76,7 +75,7 @@
!
! 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_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -91,6 +90,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_bb_mo *= 2.d0
END_PROVIDER
@ -104,7 +105,7 @@
!
! 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} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC

148
src/two_body_rdm/state_av_full_orb_2_rdm.irp.f
View File

@ -4,22 +4,16 @@
state_av_full_occ_2_rdm_ab_mo = 0.d0
integer :: i,j,k,l,iorb,jorb,korb,lorb
BEGIN_DOC
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta + beta/alpha electrons
!
! = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA
!
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
END_DOC
state_av_full_occ_2_rdm_ab_mo = 0.d0
@ -47,7 +41,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -61,7 +55,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -73,7 +67,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
enddo
enddo
@ -90,7 +84,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -104,7 +98,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -116,7 +110,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
enddo
enddo
endif
@ -135,7 +129,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -167,11 +161,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -181,8 +175,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
@ -198,11 +192,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -212,8 +206,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -231,7 +225,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -263,11 +257,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -277,8 +271,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
@ -294,11 +288,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -308,8 +302,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -328,7 +322,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -364,11 +358,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -377,8 +371,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -390,11 +384,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -403,8 +397,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -420,11 +414,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -433,8 +427,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -446,11 +440,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -459,8 +453,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -476,14 +470,14 @@
korb = list_inact(k)
! ALPHA INACTIVE - BETA ACTIVE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! BETA INACTIVE - ALPHA ACTIVE
! beta alph beta alpha
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -493,8 +487,8 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0d0
enddo
enddo
@ -510,14 +504,14 @@
korb = list_core(k)
!! BETA ACTIVE - ALPHA CORE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
!! ALPHA ACTIVE - BETA CORE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -527,8 +521,8 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0D0
enddo
enddo

2
src/two_body_rdm/test_2_rdm.irp.f
View File

@ -2,7 +2,7 @@ program test_2_rdm
implicit none
read_wf = .True.
touch read_wf
call routine_active_only
! call routine_active_only
call routine_full_mos
end

3
src/two_body_rdm/two_e_dm_mo.irp.f
View File

@ -29,7 +29,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
enddo
enddo
enddo
two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) * 2.d0
two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:)
! * 2.d0
END_PROVIDER

4
src/two_rdm_routines/davidson_like_2rdm.irp.f
View File

@ -4,8 +4,8 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
BEGIN_DOC
! if ispin == 1 :: alpha/alpha 2rdm
! == 2 :: beta /beta 2rdm
! == 3 :: alpha/beta 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
! == 3 :: alpha/beta + beta/alpha 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + ab + ba
!
! Assumes that the determinants are in psi_det
!

145
src/utils/linear_algebra.irp.f
View File

@ -1574,79 +1574,110 @@ subroutine nullify_small_elements(m,n,A,LDA,thresh)
end
subroutine restore_symmetry(m,n,A,LDA,thresh)
implicit none
BEGIN_DOC
! Tries to find the matrix elements that are the same, and sets them
! to the average value.
! If restore_symm is False, only nullify small elements
! Tries to find the matrix elements that are the same, and sets them
! to the average value.
! If restore_symm is False, only nullify small elements
END_DOC
integer, intent(in) :: m,n,LDA
double precision, intent(inout) :: A(LDA,n)
double precision, intent(in) :: thresh
integer :: i,j,k,l
logical, allocatable :: done(:,:)
double precision :: f, g, count, thresh2
double precision, allocatable :: copy(:), copy_sign(:)
integer, allocatable :: key(:), ii(:), jj(:)
integer :: sze, pi, pf, idx, i,j,k
double precision :: average, val, thresh2
thresh2 = dsqrt(thresh)
call nullify_small_elements(m,n,A,LDA,thresh)
! if (.not.restore_symm) then
! return
! endif
sze = m * n
! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)):
! - copy all values in a 1D array
! - sort 1D array
! - average nearby elements
! - for all elements, find matching value in the sorted 1D array
allocate(copy(sze),copy_sign(sze),key(sze),ii(sze),jj(sze))
allocate(done(m,n))
do j=1,n
do i=1,m
done(i,j) = A(i,j) == 0.d0
! Copy to 1D
!$OMP PARALLEL if (m>100) &
!$OMP SHARED(A,m,n,sze,copy_sign,copy,key,ii,jj) &
!$OMP PRIVATE(i,j,k) &
!$OMP DEFAULT(NONE)
!$OMP DO
do j = 1, n
do i = 1, m
k = i+(j-1)*m
copy(k) = A(i,j)
copy_sign(k) = sign(1.d0,copy(k))
copy(k) = -dabs(copy(k))
key(k) = k
ii(k) = i
jj(k) = j
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do j=1,n
do i=1,m
if ( done(i,j) ) cycle
done(i,j) = .True.
count = 1.d0
f = 1.d0/A(i,j)
do l=1,n
do k=1,m
if ( done(k,l) ) cycle
g = f * A(k,l)
if ( dabs(dabs(g) - 1.d0) < thresh2 ) then
count = count + 1.d0
if (g>0.d0) then
A(i,j) = A(i,j) + A(k,l)
else
A(i,j) = A(i,j) - A(k,l)
end if
endif
enddo
enddo
if (count > 1.d0) then
A(i,j) = A(i,j) / count
do l=1,n
do k=1,m
if ( done(k,l) ) cycle
g = f * A(k,l)
if ( dabs(dabs(g) - 1.d0) < thresh2 ) then
done(k,l) = .True.
if (g>0.d0) then
A(k,l) = A(i,j)
else
A(k,l) = -A(i,j)
end if
endif
enddo
enddo
! Sort
call dsort(copy,key,sze)
call iset_order(ii,key,sze)
call iset_order(jj,key,sze)
call dset_order(copy_sign,key,sze)
!TODO
! Parallelization with OMP
! Symmetrize
i = 1
do while (i < sze)
pi = i
pf = i
! Exit if the remaining elements are below thresh
if (-copy(i) <= thresh) exit
val = 1d0/copy(i)
do while (dabs(val * copy(pf+1) - 1d0) < thresh2)
pf = pf + 1
! if pf == sze, copy(pf+1) will not be valid
if (pf == sze) then
exit
endif
enddo
! if pi and pf are different do the average from pi to pf
if (pf - pi > 0) then
average = 0d0
do j = pi, pf
average = average + copy(j)
enddo
average = average / (pf-pi+1.d0)
do j = pi, pf
copy(j) = average
enddo
! Update i
i = pf
endif
! Update i
i = i + 1
enddo
copy(i:) = 0.d0
!$OMP PARALLEL if (sze>10000) &
!$OMP SHARED(m,sze,copy_sign,copy,key,A,ii,jj) &
!$OMP PRIVATE(i,j,k,idx) &
!$OMP DEFAULT(NONE)
! copy -> A
!$OMP DO
do k = 1, sze
i = ii(k)
j = jj(k)
A(i,j) = sign(copy(k),copy_sign(k))
enddo
!$OMP END DO
!$OMP END PARALLEL
deallocate(copy,copy_sign,key,ii,jj)
end