mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-21 11:03:29 +01:00
Merge branch 'dev-stable-tc-scf' of https://github.com/QuantumPackage/qp2 into dev-stable-tc-scf
This commit is contained in:
eginer
commit
3a6a514186
@ -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
66
codemeta.json
Normal 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
2
configure
vendored
@ -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}
|
||||
|
||||
|
||||
@ -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,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
|
||||
;;
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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)
|
||||
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -82,7 +82,7 @@
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
print*,'density and gradients provided'
|
||||
! print*,'density and gradients provided'
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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)]
|
||||
|
||||
@ -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
|
||||
|
||||
80
src/nuclei/write_pt_charges.py
Normal file
80
src/nuclei/write_pt_charges.py
Normal 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)
|
||||
|
||||
|
||||
16
src/tc_bi_ortho/print_tc_spin_dens.irp.f
Normal file
16
src/tc_bi_ortho/print_tc_spin_dens.irp.f
Normal 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
114
src/tc_bi_ortho/spin_mulliken.irp.f
Normal 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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
30
src/tc_bi_ortho/test_spin_dens.irp.f
Normal 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
|
||||
@ -183,3 +183,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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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,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
|
||||
|
||||
|
||||
@ -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,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
|
||||
!
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user