10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00

Merge pull request #264 from QuantumPackage/dev-stable

Dev stable
This commit is contained in:
Emmanuel Giner 2023-03-14 14:25:11 +01:00 committed by GitHub
commit b9b57a4f6e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
52 changed files with 17303 additions and 32 deletions

225
src/ccsd/80.ccsd_spin.bats Normal file
View File

@ -0,0 +1,225 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
source $QP_ROOT/quantum_package.rc
function run() {
thresh1=1e-6
thresh2=1e-6
test_exe scf || skip
qp set_file $1
qp edit --check
#qp run scf
qp set_frozen_core
qp set utils_cc cc_par_t true
qp set utils_cc cc_thresh_conv 1e-12
file="$(echo $1 | sed 's/.ezfio//g')"
qp run ccsd_spin_orb | tee $file.ccsd.out
energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
#rm $file.ccsd.out
eq $energy1 $2 $thresh1
eq $energy2 $3 $thresh2
}
@test "b2_stretched" {
run b2_stretched.ezfio -49.136487344382 -0.003497589175
}
@test "be" {
run be.ezfio -14.623559003577 -0.000230982022
}
@test "c2h2" {
run c2h2.ezfio -12.394008897618 -0.010790491561
}
@test "ch4" {
run ch4.ezfio -40.390721785799 -0.004476100282
}
@test "clf" {
run clf.ezfio -559.186562904081 -0.006577143392
}
@test "clo" {
run clo.ezfio -534.564874409332 -0.007584571424
}
@test "co2" {
run co2.ezfio -188.129602527766 -0.018040668885
}
@test "dhno" {
run dhno.ezfio -130.816650109473 -0.012197331453
}
@test "f2" {
run f2.ezfio -199.287826338097 -0.017592872692
}
@test "f" {
run f.ezfio -99.616644511121 -0.003624525307
}
@test "h2o2" {
run h2o2.ezfio -151.182552729963 -0.009511682086
}
@test "h2o" {
run h2o.ezfio -76.237710276526 -0.003001800577
}
@test "h2s" {
run h2s.ezfio -398.861214015390 -0.003300559757
}
@test "h3coh" {
run h3coh.ezfio -115.221296424969 -0.003566171432
}
@test "hbo" {
run hbo.ezfio -100.213539770415 -0.006851489212
}
@test "hcn" {
run hcn.ezfio -93.190247992657 -0.013418135043
}
@test "hco" {
run hco.ezfio -113.405413962350 -0.007973455337
}
@test "lif" {
run lif.ezfio -107.270402903250 -0.007742969005
}
@test "n2" {
run n2.ezfio -109.355358930472 -0.018477744342
}
@test "n2h4" {
run n2h4.ezfio -111.556885923139 -0.009048077008
}
@test "nh3" {
run nh3.ezfio -56.465503060954 -0.007638273755
}
@test "oh" {
run oh.ezfio -75.614606132774 -0.004126661739
}
@test "sih2_3b1" {
run sih2_3b1.ezfio -290.016780973072 -0.000497825874
}
@test "sih3" {
run sih3.ezfio -5.575343504534 -0.002094123268
}
@test "so" {
run so.ezfio -26.035945178665 -0.010594351274
}
#@test "b2_stretched" {
#run b2_stretched.ezfio -49.136487344382 -49.139984933557
#}
#
#@test "be" {
#run be.ezfio -14.623559003577 -14.623789985599
#}
#
#@test "c2h2" {
#run c2h2.ezfio -12.394008897618 -12.404799389179
#}
#
#@test "ch4" {
#run ch4.ezfio -40.390721784961 -40.395197884406
#}
#
#@test "clf" {
#run clf.ezfio -559.186562906072 -559.193140046904
#}
#
#@test "clo" {
#run clo.ezfio -534.564874409333 -534.572458980757
#}
#
#@test "co2" {
#run co2.ezfio -188.129602511724 -188.147643198675
#}
#
#@test "dhno" {
#run dhno.ezfio -130.816650109473 -130.828847440925
#}
#
#@test "f2" {
#run f2.ezfio -199.287826338097 -199.305419210789
#}
#
#@test "f" {
#run f.ezfio -99.616644511120 -99.620269036428
#}
#
#@test "h2o2" {
#run h2o2.ezfio -151.182552729963 -151.192064412049
#}
#
#@test "h2o" {
#run h2o.ezfio -76.237710276526 -76.240712077103
#}
#
#@test "h2s" {
#run h2s.ezfio -398.861214015416 -398.864514575146
#}
#
#@test "h3coh" {
#run h3coh.ezfio -115.221296424969 -115.224862596401
#}
#
#@test "hbo" {
#run hbo.ezfio -100.213539770415 -100.220391259627
#}
#
#@test "hcn" {
#run hcn.ezfio -93.190247983000 -93.203666131216
#}
#
#@test "hco" {
#run hco.ezfio -113.405413962350 -113.413387417687
#}
#
#@test "lif" {
#run lif.ezfio -107.270402903211 -107.278145872216
#}
#
#@test "n2" {
#run n2.ezfio -109.355358930472 -109.373836674814
#}
#
#@test "n2h4" {
#run n2h4.ezfio -111.556885922642 -111.565934000556
#}
#
#@test "nh3" {
#run nh3.ezfio -56.465503060954 -56.473141334709
#}
#
#@test "oh" {
#run oh.ezfio -75.614606131897 -75.618732794235
#}
#
#@test "sih2_3b1" {
#run sih2_3b1.ezfio -290.016780973071 -290.017278798946
#}
#
#@test "sih3" {
#run sih3.ezfio -5.575343504534 -5.577437627802
#}
#
#@test "so" {
#run so.ezfio -26.035945181998 -26.046539528491
#}

225
src/ccsd/81.ccsd_space.bats Normal file
View File

@ -0,0 +1,225 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
source $QP_ROOT/quantum_package.rc
function run() {
thresh1=1e-6
thresh2=1e-6
test_exe scf || skip
qp set_file $1
qp edit --check
#qp run scf
qp set_frozen_core
qp set utils_cc cc_par_t true
qp set utils_cc cc_thresh_conv 1e-12
file="$(echo $1 | sed 's/.ezfio//g')"
qp run ccsd_space_orb | tee $file.ccsd.out
energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
#rm $file.ccsd.out
eq $energy1 $2 $thresh1
eq $energy2 $3 $thresh2
}
@test "b2_stretched" {
run b2_stretched.ezfio -49.136487344382 -0.003497589175
}
@test "be" {
run be.ezfio -14.623559003577 -0.000230982022
}
@test "c2h2" {
run c2h2.ezfio -12.394008897618 -0.010790491561
}
@test "ch4" {
run ch4.ezfio -40.390721785799 -0.004476100282
}
@test "clf" {
run clf.ezfio -559.186562904081 -0.006577143392
}
#@test "clo" {
#run clo.ezfio -534.564874409332 -0.007584571424
#}
@test "co2" {
run co2.ezfio -188.129602527766 -0.018040668885
}
#@test "dhno" {
#run dhno.ezfio -130.816650109473 -0.012197331453
#}
@test "f2" {
run f2.ezfio -199.287826338097 -0.017592872692
}
#@test "f" {
#run f.ezfio -99.616644511121 -0.003624525307
#}
@test "h2o2" {
run h2o2.ezfio -151.182552729963 -0.009511682086
}
@test "h2o" {
run h2o.ezfio -76.237710276526 -0.003001800577
}
@test "h2s" {
run h2s.ezfio -398.861214015390 -0.003300559757
}
@test "h3coh" {
run h3coh.ezfio -115.221296424969 -0.003566171432
}
@test "hbo" {
run hbo.ezfio -100.213539770415 -0.006851489212
}
@test "hcn" {
run hcn.ezfio -93.190247992657 -0.013418135043
}
#@test "hco" {
#run hco.ezfio -113.405413962350 -0.007973455337
#}
@test "lif" {
run lif.ezfio -107.270402903250 -0.007742969005
}
@test "n2" {
run n2.ezfio -109.355358930472 -0.018477744342
}
@test "n2h4" {
run n2h4.ezfio -111.556885923139 -0.009048077008
}
@test "nh3" {
run nh3.ezfio -56.465503060954 -0.007638273755
}
#@test "oh" {
#run oh.ezfio -75.614606132774 -0.004126661739
#}
#@test "sih2_3b1" {
#run sih2_3b1.ezfio -290.016780973072 -0.000497825874
#}
#@test "sih3" {
#run sih3.ezfio -5.575343504534 -0.002094123268
#}
#@test "so" {
#run so.ezfio -26.035945178665 -0.010594351274
#}
#@test "b2_stretched" {
#run b2_stretched.ezfio -49.136487344382 -49.139984933557
#}
#
#@test "be" {
#run be.ezfio -14.623559003577 -14.623789985599
#}
#
#@test "c2h2" {
#run c2h2.ezfio -12.394008897618 -12.404799389179
#}
#
#@test "ch4" {
#run ch4.ezfio -40.390721784961 -40.395197884406
#}
#
#@test "clf" {
#run clf.ezfio -559.186562906072 -559.193140046904
#}
#
##@test "clo" {
##run clo.ezfio -534.564874409333 -534.572458980757
##}
#
#@test "co2" {
#run co2.ezfio -188.129602511724 -188.147643198675
#}
#
##@test "dhno" {
##run dhno.ezfio -130.816650109473 -130.828847440925
##}
#
#@test "f2" {
#run f2.ezfio -199.287826338097 -199.305419210789
#}
#
##@test "f" {
##run f.ezfio -99.616644511120 -99.620269036428
##}
#
#@test "h2o2" {
#run h2o2.ezfio -151.182552729963 -151.192064412049
#}
#
#@test "h2o" {
#run h2o.ezfio -76.237710276526 -76.240712077103
#}
#
#@test "h2s" {
#run h2s.ezfio -398.861214015416 -398.864514575146
#}
#
#@test "h3coh" {
#run h3coh.ezfio -115.221296424969 -115.224862596401
#}
#
#@test "hbo" {
#run hbo.ezfio -100.213539770415 -100.220391259627
#}
#
#@test "hcn" {
#run hcn.ezfio -93.190247983000 -93.203666131216
#}
#
##@test "hco" {
##run hco.ezfio -113.405413962350 -113.413387417687
##}
#
#@test "lif" {
#run lif.ezfio -107.270402903211 -107.278145872216
#}
#
#@test "n2" {
#run n2.ezfio -109.355358930472 -109.373836674814
#}
#
#@test "n2h4" {
#run n2h4.ezfio -111.556885922642 -111.565934000556
#}
#
#@test "nh3" {
#run nh3.ezfio -56.465503060954 -56.473141334709
#}
#
##@test "oh" {
##run oh.ezfio -75.614606131897 -75.618732794235
##}
#
##@test "sih2_3b1" {
##run sih2_3b1.ezfio -290.016780973071 -290.017278798946
##}
#
##@test "sih3" {
##run sih3.ezfio -5.575343504534 -5.577437627802
##}
#
##@test "so" {
##run so.ezfio -26.035945181998 -26.046539528491
##}

2
src/ccsd/NEED Normal file
View File

@ -0,0 +1,2 @@
hartree_fock
utils_cc

31
src/ccsd/README.md Normal file
View File

@ -0,0 +1,31 @@
# CCSD in spin orbitals and spatial orbitals
CCSD and CCSD(T) in spin orbitals for open and closed shell systems.
CCSD and CCSD(T) in spatial orbitals for closed shell systems.
## Calculations
The program will automatically choose the version in spin or spatial orbitals
To run the general program:
```
qp run ccsd
```
Nevertheless, you can enforce the run in spin orbitals with
```
qp run ccsd_spin_orb
```
## Settings
The settings can be changed with:
```
qp set utils_cc cc_#param #val
```
For more informations on the settings, look at the module utils_cc and its documentation.
## Org files
The org files are stored in the directory org in order to avoid overwriting on user changes.
The org files can be modified, to export the change to the source code, run
```
./TANGLE_org_mode.sh and
mv *.irp.f ../.
```

18
src/ccsd/ccsd.irp.f Normal file
View File

@ -0,0 +1,18 @@
program ccsd
implicit none
BEGIN_DOC
! CCSD program
END_DOC
read_wf = .True.
touch read_wf
if (.not. cc_ref_is_open_shell) then
call run_ccsd_space_orb
else
call run_ccsd_spin_orb
endif
end

View File

@ -0,0 +1,12 @@
! Code
program ccsd
implicit none
read_wf = .True.
touch read_wf
call run_ccsd_space_orb
end

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,16 @@
! Prog
program ccsd
implicit none
BEGIN_DOC
! CCSD in spin orbitals
END_DOC
read_wf = .True.
touch read_wf
call run_ccsd_spin_orb
end

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,412 @@
! Dumb way
subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: energy
double precision, allocatable :: W(:,:,:,:,:,:)
double precision, allocatable :: V(:,:,:,:,:,:)
integer :: i,j,k,a,b,c
allocate(W(nO,nO,nO,nV,nV,nV))
allocate(V(nO,nO,nO,nV,nV,nV))
call form_w(nO,nV,t2,W)
call form_v(nO,nV,t1,W,V)
energy = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c)) !delta_ooovvv(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
energy = energy / 3d0
deallocate(V,W)
end
subroutine form_w(nO,nV,t2,W)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
integer :: i,j,k,l,a,b,c,d
W = 0d0
do c = 1, nV
print*,'W:',c,'/',nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do d = 1, nV
W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
! chem (bd|ai)
! phys <ba|di>
+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
enddo
do l = 1, nO
W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
! chem (ck|jl)
! phys <cj|kl>
- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
enddo
enddo
enddo
enddo
enddo
enddo
enddo
end
subroutine form_v(nO,nV,t1,w,v)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: W(nO, nO, nO, nV, nV, nV)
double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
integer :: i,j,k,a,b,c
V = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
enddo
enddo
enddo
enddo
enddo
enddo
end
! Main
subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
double precision, intent(out) :: energy
double precision, allocatable :: W(:,:,:,:,:,:)
double precision, allocatable :: V(:,:,:,:,:,:)
double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
integer :: i,j,k,l,a,b,c,d
double precision :: e,ta,tb, delta, delta_ijk
!allocate(W(nV,nV,nV,nO,nO,nO))
!allocate(V(nV,nV,nV,nO,nO,nO))
allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
! Temporary arrays
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
!$OMP PRIVATE(a,b,c,d,i,j,k,l) &
!$OMP DEFAULT(NONE)
!v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
!$OMP DO collapse(3)
do i = 1, nO
do a = 1, nV
do b = 1, nV
do d = 1, nV
X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
do j = 1, nO
do k = 1, nO
do c = 1, nV
do d = 1, nV
T_vvoo(d,c,k,j) = t2(k,j,c,d)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vooo(c,j,k,l) * t2(i,l,a,b) &
!X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
!$OMP DO collapse(3)
do k = 1, nO
do j = 1, nO
do c = 1, nV
do l = 1, nO
X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
do i = 1, nO
do b = 1, nV
do a = 1, nV
do l = 1, nO
T_ovvo(l,a,b,i) = t2(i,l,a,b)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vvoo(b,c,j,k) * t1(i,a) &
!X_vvoo(b,c,k,j) * T1_vo(a,i) &
!$OMP DO collapse(3)
do j = 1, nO
do k = 1, nO
do c = 1, nV
do b = 1, nV
X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(1)
do i = 1, nO
do a = 1, nV
T_vo(a,i) = t1(i,a)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(ta)
energy = 0d0
do i = 1, nO
do j = 1, nO
do k = 1, nO
delta_ijk = f_o(i) + f_o(j) + f_o(k)
call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
!$OMP PARALLEL &
!$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
!$OMP PRIVATE(a,b,c,e,delta) &
!$OMP DEFAULT(NONE)
e = 0d0
!$OMP DO
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
!energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c)) !delta_ooovvv(i,j,k,a,b,c)
e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
* (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
enddo
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
energy = energy + e
!$OMP END CRITICAL
!$OMP END PARALLEL
enddo
enddo
call wall_time(tb)
write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
enddo
energy = energy / 3d0
deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
!deallocate(V,W)
end
! W_ijk
subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
implicit none
integer, intent(in) :: nO,nV,i,j,k
!double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
double precision, intent(in) :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
integer :: l,a,b,c,d
!W = 0d0
!do i = 1, nO
! do j = 1, nO
! do k = 1, nO
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
!$OMP PRIVATE(a,b,c,d,l) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
W(a,b,c) = 0d0
do d = 1, nV
!W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
W(a,b,c) = W(a,b,c) &
! chem (bd|ai)
! phys <ba|di>
!+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
!+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
!+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
!+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
!+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+ X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+ X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+ X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+ X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+ X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+ X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
do l = 1, nO
!W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
W(a,b,c) = W(a,b,c) &
! chem (ck|jl)
! phys <cj|kl>
!- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
!- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
!- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
!- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
!- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
!- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
- X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
- X_ovoo(l,b,k,j) * T_ovvo(l,a,c,i) & ! bc kj
- X_ovoo(l,b,i,j) * T_ovvo(l,c,a,k) & ! prev ac ik
- X_ovoo(l,a,j,i) * T_ovvo(l,c,b,k) & ! prev ab ij
- X_ovoo(l,a,k,i) * T_ovvo(l,b,c,j) & ! prev bc kj
- X_ovoo(l,c,i,k) * T_ovvo(l,b,a,j) ! prev ac ik
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! enddo
! enddo
!enddo
end
! V_ijk
subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
implicit none
integer, intent(in) :: nO,nV,i,j,k
!double precision, intent(in) :: t1(nO,nV)
double precision, intent(in) :: T_vo(nV,nO)
double precision, intent(in) :: X_vvoo(nV,nV,nO,nO)
double precision, intent(in) :: W(nV,nV,nV)!,nO,nO,nO)
double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
integer :: a,b,c
!V = 0d0
!do i = 1, nO
! do j = 1, nO
! do k = 1, nO
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
!$OMP PRIVATE(a,b,c) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
!V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
V(a,b,c) = W(a,b,c) &
!+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
!+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
!+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+ X_vvoo(b,c,k,j) * T_vo(a,i) &
+ X_vvoo(a,c,k,i) * T_vo(b,j) &
+ X_vvoo(a,b,j,i) * T_vo(c,k)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! enddo
! enddo
!enddo
end

View File

@ -0,0 +1,376 @@
! v1
subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
implicit none
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
double precision, intent(in) :: v_ooov(nO,nO,nO,nV)
double precision, intent(in) :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
double precision, intent(out) :: energy
double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
double precision :: e_t, e_st, e_dt, delta_abc, delta
integer :: i,j,k,l,m,a,b,c,d,e
allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
t3 = 0d0
! T3
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
do e = 1, nV
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+ t2(j,i,c,e) * v_vvvo(b,a,e,k) ! + P(ik) P(ac)
enddo
do m = 1, nO
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+ t2(m,k,b,a) * v_ooov(j,i,m,c) ! + P(ik) P(ac)
enddo
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
enddo
enddo
enddo
enddo
enddo
enddo
! E_T
e_t = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_t = e_t / 36d0
! E_ST
s = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_st = 0d0
do a = 1, nV
do i = 1, nO
e_st = e_st + s(i,a) * t1(i,a)
enddo
enddo
e_st = e_st * 0.25d0
! E_DT
e_dt = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_dt = e_dt * 0.25d0
! (T)
!print*,e_t,e_st,e_dt
energy = e_t + e_st + e_dt
deallocate(t3,s)
end
! v2
subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
implicit none
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
double precision, intent(in) :: v_ooov(nO,nO,nO,nV)
double precision, intent(in) :: v_vvoo(nV,nV,nO,nO)
double precision, intent(out) :: energy
double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
double precision :: e_st, delta_abc, delta, ta, tb
integer :: i,j,k,l,m,a,b,c,d,e
allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
nV,nV,nV,nO, v_vvvo)
! Init
s = 0d0
e_t = 0d0
e_st = 0d0
e_dt = 0d0
call wall_time(ta)
!$OMP PARALLEL &
!$OMP PRIVATE(i,j,k,m,a,b,c,e) &
!$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
!$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do c = 1, nV
do b = 1, nV
do i = 1, nO
do e = 1, nV
A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$omp do collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do e = 1, nV
T_voov(e,j,k,a) = t2(j,k,a,e)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do m = 1, nO
B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
enddo
enddo
enddo
enddo
!$omp end do
do c = 1, nV
do b = 1, nV
! T3(:,:,:,:,b,c)
! Init
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
t3_bc(i,j,k,a) = 0d0
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do e = 1, nV
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
!+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
!- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
!- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
!- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
!- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
!+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
!+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
!+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
!+ t2(j,i,c,e) * v_vvvo(b,a,e,k) ! + P(ik) P(ac)
+ T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
- T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
- T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
- T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
- T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+ T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+ T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+ T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+ T_voov(e,j,i,c) * A_vovv(e,k,b,a) ! + P(ik) P(ac)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do m = 1, nO
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
!+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
!- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
!- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
!- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
!- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
!+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
!+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
!+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
!+ t2(m,k,b,a) * v_ooov(j,i,m,c) ! + P(ik) P(ac)
+ t2(m,i,b,c) * B_ooov(m,j,k,a) &
- t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
- t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
- t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
- t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+ t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+ t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+ t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+ t2(m,k,b,a) * B_ooov(m,j,i,c) ! + P(ik) P(ac)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
enddo
enddo
enddo
enddo
!$OMP END DO
! E_T
!$OMP DO
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
! E_ST
!$OMP DO
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
! E_DT
!$OMP DO
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO
enddo
!$OMP MASTER
call wall_time(tb)
write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
!$OMP END MASTER
enddo
!$OMP END PARALLEL
do a = 2, nV
e_t(1) = e_t(1) + e_t(a)
enddo
do a = 2, nV
e_dt(1) = e_dt(1) + e_dt(a)
enddo
e_t = e_t / 36d0
do a = 1, nV
do i = 1, nO
e_st = e_st + s(i,a) * t1(i,a)
enddo
enddo
e_st = e_st * 0.25d0
e_dt = e_dt * 0.25d0
! (T)
!print*,e_t(1),e_st,e_dt(1)
energy = e_t(1) + e_st + e_dt(1)
deallocate(t3_bc,s)
end

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,428 @@
Ref:
Integral-Direct and Parallel Implementation of the CCSD(T) Method:
Algorithmic Developments and Large-Scale Applications
László Gyevi-Nagy, Mihály Kállay, and Péter R. Nagy
J. Chem. Theory Comput. 2020, 16, 1, 366384
https://doi.org/10.1021/acs.jctc.9b00957
* Dumb way
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: energy
double precision, allocatable :: W(:,:,:,:,:,:)
double precision, allocatable :: V(:,:,:,:,:,:)
integer :: i,j,k,a,b,c
allocate(W(nO,nO,nO,nV,nV,nV))
allocate(V(nO,nO,nO,nV,nV,nV))
call form_w(nO,nV,t2,W)
call form_v(nO,nV,t1,W,V)
energy = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c)) !delta_ooovvv(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
energy = energy / 3d0
deallocate(V,W)
end
#+END_SRC
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine form_w(nO,nV,t2,W)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
integer :: i,j,k,l,a,b,c,d
W = 0d0
do c = 1, nV
print*,'W:',c,'/',nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do d = 1, nV
W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
! chem (bd|ai)
! phys <ba|di>
+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
enddo
do l = 1, nO
W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
! chem (ck|jl)
! phys <cj|kl>
- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
enddo
enddo
enddo
enddo
enddo
enddo
enddo
end
#+END_SRC
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine form_v(nO,nV,t1,w,v)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: W(nO, nO, nO, nV, nV, nV)
double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
integer :: i,j,k,a,b,c
V = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
enddo
enddo
enddo
enddo
enddo
enddo
end
#+END_SRC
* Better way
** Main
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
double precision, intent(out) :: energy
double precision, allocatable :: W(:,:,:,:,:,:)
double precision, allocatable :: V(:,:,:,:,:,:)
double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
integer :: i,j,k,l,a,b,c,d
double precision :: e,ta,tb, delta, delta_ijk
!allocate(W(nV,nV,nV,nO,nO,nO))
!allocate(V(nV,nV,nV,nO,nO,nO))
allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
! Temporary arrays
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
!$OMP PRIVATE(a,b,c,d,i,j,k,l) &
!$OMP DEFAULT(NONE)
!v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
!$OMP DO collapse(3)
do i = 1, nO
do a = 1, nV
do b = 1, nV
do d = 1, nV
X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
do j = 1, nO
do k = 1, nO
do c = 1, nV
do d = 1, nV
T_vvoo(d,c,k,j) = t2(k,j,c,d)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vooo(c,j,k,l) * t2(i,l,a,b) &
!X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
!$OMP DO collapse(3)
do k = 1, nO
do j = 1, nO
do c = 1, nV
do l = 1, nO
X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
do i = 1, nO
do b = 1, nV
do a = 1, nV
do l = 1, nO
T_ovvo(l,a,b,i) = t2(i,l,a,b)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vvoo(b,c,j,k) * t1(i,a) &
!X_vvoo(b,c,k,j) * T1_vo(a,i) &
!$OMP DO collapse(3)
do j = 1, nO
do k = 1, nO
do c = 1, nV
do b = 1, nV
X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(1)
do i = 1, nO
do a = 1, nV
T_vo(a,i) = t1(i,a)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(ta)
energy = 0d0
do i = 1, nO
do j = 1, nO
do k = 1, nO
delta_ijk = f_o(i) + f_o(j) + f_o(k)
call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
!$OMP PARALLEL &
!$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
!$OMP PRIVATE(a,b,c,e,delta) &
!$OMP DEFAULT(NONE)
e = 0d0
!$OMP DO
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
!energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c)) !delta_ooovvv(i,j,k,a,b,c)
e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
* (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
enddo
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
energy = energy + e
!$OMP END CRITICAL
!$OMP END PARALLEL
enddo
enddo
call wall_time(tb)
write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
enddo
energy = energy / 3d0
deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
!deallocate(V,W)
end
#+END_SRC
** W_ijk
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
implicit none
integer, intent(in) :: nO,nV,i,j,k
!double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
double precision, intent(in) :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
integer :: l,a,b,c,d
!W = 0d0
!do i = 1, nO
! do j = 1, nO
! do k = 1, nO
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
!$OMP PRIVATE(a,b,c,d,l) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
W(a,b,c) = 0d0
do d = 1, nV
!W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
W(a,b,c) = W(a,b,c) &
! chem (bd|ai)
! phys <ba|di>
!+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
!+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
!+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
!+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
!+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+ X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+ X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+ X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+ X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+ X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+ X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
do l = 1, nO
!W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
W(a,b,c) = W(a,b,c) &
! chem (ck|jl)
! phys <cj|kl>
!- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
!- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
!- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
!- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
!- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
!- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
- X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
- X_ovoo(l,b,k,j) * T_ovvo(l,a,c,i) & ! bc kj
- X_ovoo(l,b,i,j) * T_ovvo(l,c,a,k) & ! prev ac ik
- X_ovoo(l,a,j,i) * T_ovvo(l,c,b,k) & ! prev ab ij
- X_ovoo(l,a,k,i) * T_ovvo(l,b,c,j) & ! prev bc kj
- X_ovoo(l,c,i,k) * T_ovvo(l,b,a,j) ! prev ac ik
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! enddo
! enddo
!enddo
end
#+END_SRC
** V_ijk
#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
implicit none
integer, intent(in) :: nO,nV,i,j,k
!double precision, intent(in) :: t1(nO,nV)
double precision, intent(in) :: T_vo(nV,nO)
double precision, intent(in) :: X_vvoo(nV,nV,nO,nO)
double precision, intent(in) :: W(nV,nV,nV)!,nO,nO,nO)
double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
integer :: a,b,c
!V = 0d0
!do i = 1, nO
! do j = 1, nO
! do k = 1, nO
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
!$OMP PRIVATE(a,b,c) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(2)
do c = 1, nV
do b = 1, nV
do a = 1, nV
!V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
V(a,b,c) = W(a,b,c) &
!+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
!+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
!+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+ X_vvoo(b,c,k,j) * T_vo(a,i) &
+ X_vvoo(a,c,k,i) * T_vo(b,j) &
+ X_vvoo(a,b,j,i) * T_vo(c,k)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! enddo
! enddo
!enddo
end
#+END_SRC

View File

@ -0,0 +1,385 @@
* CCSD(T) spin orb
Ref:
John D. Watts, Jürgen Gauss, and Rodney J. Bartlett
J. Chem. Phys. 98, 8718 (1993)
http://dx.doi.org/10.1063/1.464480
** v1
#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
implicit none
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
double precision, intent(in) :: v_ooov(nO,nO,nO,nV)
double precision, intent(in) :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
double precision, intent(out) :: energy
double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
double precision :: e_t, e_st, e_dt, delta_abc, delta
integer :: i,j,k,l,m,a,b,c,d,e
allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
t3 = 0d0
! T3
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
do e = 1, nV
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+ t2(j,i,c,e) * v_vvvo(b,a,e,k) ! + P(ik) P(ac)
enddo
do m = 1, nO
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+ t2(m,k,b,a) * v_ooov(j,i,m,c) ! + P(ik) P(ac)
enddo
t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
enddo
enddo
enddo
enddo
enddo
enddo
! E_T
e_t = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_t = e_t / 36d0
! E_ST
s = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_st = 0d0
do a = 1, nV
do i = 1, nO
e_st = e_st + s(i,a) * t1(i,a)
enddo
enddo
e_st = e_st * 0.25d0
! E_DT
e_dt = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
enddo
enddo
enddo
enddo
enddo
enddo
e_dt = e_dt * 0.25d0
! (T)
!print*,e_t,e_st,e_dt
energy = e_t + e_st + e_dt
deallocate(t3,s)
end
#+end_src
** v2
#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
implicit none
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
double precision, intent(in) :: v_ooov(nO,nO,nO,nV)
double precision, intent(in) :: v_vvoo(nV,nV,nO,nO)
double precision, intent(out) :: energy
double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
double precision :: e_st, delta_abc, delta, ta, tb
integer :: i,j,k,l,m,a,b,c,d,e
allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
nV,nV,nV,nO, v_vvvo)
! Init
s = 0d0
e_t = 0d0
e_st = 0d0
e_dt = 0d0
call wall_time(ta)
!$OMP PARALLEL &
!$OMP PRIVATE(i,j,k,m,a,b,c,e) &
!$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
!$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do c = 1, nV
do b = 1, nV
do i = 1, nO
do e = 1, nV
A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$omp do collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do e = 1, nV
T_voov(e,j,k,a) = t2(j,k,a,e)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do m = 1, nO
B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
enddo
enddo
enddo
enddo
!$omp end do
do c = 1, nV
do b = 1, nV
! T3(:,:,:,:,b,c)
! Init
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
t3_bc(i,j,k,a) = 0d0
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do e = 1, nV
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
!+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
!- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
!- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
!- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
!- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
!+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
!+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
!+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
!+ t2(j,i,c,e) * v_vvvo(b,a,e,k) ! + P(ik) P(ac)
+ T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
- T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
- T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
- T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
- T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+ T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+ T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+ T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+ T_voov(e,j,i,c) * A_vovv(e,k,b,a) ! + P(ik) P(ac)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO collapse(3)
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
do m = 1, nO
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
!+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
!- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
!- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
!- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
!- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
!+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
!+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
!+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
!+ t2(m,k,b,a) * v_ooov(j,i,m,c) ! + P(ik) P(ac)
+ t2(m,i,b,c) * B_ooov(m,j,k,a) &
- t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
- t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
- t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
- t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+ t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+ t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+ t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+ t2(m,k,b,a) * B_ooov(m,j,i,c) ! + P(ik) P(ac)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
enddo
enddo
enddo
enddo
!$OMP END DO
! E_T
!$OMP DO
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
! E_ST
!$OMP DO
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
! E_DT
!$OMP DO
do a = 1, nV
do k = 1, nO
do j = 1, nO
do i = 1, nO
e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
enddo
enddo
enddo
enddo
!$OMP END DO
enddo
!$OMP MASTER
call wall_time(tb)
write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
!$OMP END MASTER
enddo
!$OMP END PARALLEL
do a = 2, nV
e_t(1) = e_t(1) + e_t(a)
enddo
do a = 2, nV
e_dt(1) = e_dt(1) + e_dt(a)
enddo
e_t = e_t / 36d0
do a = 1, nV
do i = 1, nO
e_st = e_st + s(i,a) * t1(i,a)
enddo
enddo
e_st = e_st * 0.25d0
e_dt = e_dt * 0.25d0
! (T)
!print*,e_t(1),e_st,e_dt(1)
energy = e_t(1) + e_st + e_dt(1)
deallocate(t3_bc,s)
end
#+end_src

View File

@ -3,3 +3,4 @@ zmq
mpi mpi
iterations iterations
csf csf
mol_properties

View File

@ -108,6 +108,7 @@ subroutine run_cipsi
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
call print_mol_properties()
N_iter += 1 N_iter += 1
if (qp_stop()) exit if (qp_stop()) exit
@ -156,6 +157,7 @@ subroutine run_cipsi
pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
call print_mol_properties()
endif endif
call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err) call pt2_dealloc(pt2_data_err)

View File

@ -98,6 +98,7 @@ subroutine run_stochastic_cipsi
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
call print_mol_properties()
N_iter += 1 N_iter += 1
if (qp_stop()) exit if (qp_stop()) exit
@ -136,6 +137,7 @@ subroutine run_stochastic_cipsi
pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2) pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
call print_mol_properties()
endif endif
call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err) call pt2_dealloc(pt2_data_err)

View File

@ -10,6 +10,7 @@ function run() {
qp set determinants n_states 2 qp set determinants n_states 2
qp set davidson threshold_davidson 1.e-12 qp set davidson threshold_davidson 1.e-12
qp set davidson n_states_diag 24 qp set davidson n_states_diag 24
qp run cis
qp run cisd qp run cisd
energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)" energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)" energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
@ -20,26 +21,31 @@ function run() {
@test "B-B" { # @test "B-B" { #
qp set_file b2_stretched.ezfio qp set_file b2_stretched.ezfio
qp set_frozen_core
run -49.120607088648597 -49.055152453388231 run -49.120607088648597 -49.055152453388231
} }
@test "SiH2_3B1" { # 1.53842s 3.53856s @test "SiH2_3B1" { # 1.53842s 3.53856s
qp set_file sih2_3b1.ezfio qp set_file sih2_3b1.ezfio
qp set_frozen_core
run -290.015949171697 -289.805036176618 run -290.015949171697 -289.805036176618
} }
@test "HBO" { # 4.42968s 19.6099s @test "HBO" { # 4.42968s 19.6099s
qp set_file hbo.ezfio qp set_file hbo.ezfio
qp set_frozen_core
run -100.2019254455993 -99.79484127741013 run -100.2019254455993 -99.79484127741013
} }
@test "HCO" { # 6.6077s 28.6801s @test "HCO" { # 6.6077s 28.6801s
qp set_file hco.ezfio qp set_file hco.ezfio
qp set_frozen_core
run -113.39088802205114 -113.22204293108558 run -113.39088802205114 -113.22204293108558
} }
@test "H2O" { # 7.0651s 30.6642s @test "H2O" { # 7.0651s 30.6642s
qp set_file h2o.ezfio qp set_file h2o.ezfio
qp set_frozen_core
run -76.22975602077072 -75.80609108747208 run -76.22975602077072 -75.80609108747208
} }
@ -50,6 +56,7 @@ function run() {
@test "H2S" { # 7.42152s 32.5461s @test "H2S" { # 7.42152s 32.5461s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file h2s.ezfio qp set_file h2s.ezfio
qp set_frozen_core
run -398.853701416768 -398.519020035337 run -398.853701416768 -398.519020035337
} }
@ -70,6 +77,7 @@ function run() {
@test "OH" { # 18.2159s 1.28453m @test "OH" { # 18.2159s 1.28453m
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file oh.ezfio qp set_file oh.ezfio
qp set_frozen_core
run -75.6087472926588 -75.5370393736601 run -75.6087472926588 -75.5370393736601
} }
@ -83,6 +91,7 @@ function run() {
@test "SiH3" { # 20.2202s 1.38648m @test "SiH3" { # 20.2202s 1.38648m
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file sih3.ezfio qp set_file sih3.ezfio
qp set_frozen_core
run -5.57096611856522 -5.30950347928823 run -5.57096611856522 -5.30950347928823
} }
@ -103,6 +112,7 @@ function run() {
@test "H3COH" { # 24.7248s 1.85043m @test "H3COH" { # 24.7248s 1.85043m
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file h3coh.ezfio qp set_file h3coh.ezfio
qp set_frozen_core
run -115.204958752377 -114.755913828245 run -115.204958752377 -114.755913828245
} }
@ -117,6 +127,7 @@ function run() {
@test "ClF" { # 30.3225s @test "ClF" { # 30.3225s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file clf.ezfio qp set_file clf.ezfio
qp set_frozen_core
run -559.162476603880 -558.792395927088 run -559.162476603880 -558.792395927088
} }
@ -130,6 +141,7 @@ function run() {
@test "ClO" { # 37.6949s @test "ClO" { # 37.6949s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file clo.ezfio qp set_file clo.ezfio
qp set_frozen_core
run -534.5404021326773 -534.3818725793897 run -534.5404021326773 -534.3818725793897
} }
@ -150,6 +162,7 @@ function run() {
@test "SO" { # 51.2476s @test "SO" { # 51.2476s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file so.ezfio qp set_file so.ezfio
qp set_frozen_core
run -26.0131812819785 -25.7053111980226 run -26.0131812819785 -25.7053111980226
} }

View File

@ -69,7 +69,9 @@ subroutine run
do i = 1,N_states do i = 1,N_states
k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1) k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
delta_E = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int) delta_E = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2) if (elec_alpha_num + elec_beta_num >= 4) then
cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
endif
enddo enddo
print *, 'N_det = ', N_det print *, 'N_det = ', N_det
print*,'' print*,''
@ -78,26 +80,43 @@ subroutine run
do i = 1,N_states do i = 1,N_states
print *, i, CI_energy(i) print *, i, CI_energy(i)
enddo enddo
print*,'' if (elec_alpha_num + elec_beta_num >= 4) then
print*,'******************************' print*,''
print *, 'CISD+Q Energies' print*,'******************************'
do i = 1,N_states print *, 'CISD+Q Energies'
print *, i, cisdq(i) do i = 1,N_states
enddo print *, i, cisdq(i)
enddo
endif
if (N_states > 1) then if (N_states > 1) then
print*,'' if (elec_alpha_num + elec_beta_num >= 4) then
print*,'******************************' print*,''
print*,'Excitation energies (au) (CISD+Q)' print*,'******************************'
do i = 2, N_states print*,'Excitation energies (au) (CISD+Q)'
print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1) do i = 2, N_states
enddo print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1)
print*,'' enddo
print*,'******************************' print*,''
print*,'Excitation energies (eV) (CISD+Q)' print*,'******************************'
do i = 2, N_states print*,'Excitation energies (eV) (CISD+Q)'
print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, & do i = 2, N_states
(cisdq(i) - cisdq(1)) / 0.0367502d0 print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev, &
enddo (cisdq(i) - cisdq(1)) * ha_to_ev
enddo
else
print*,''
print*,'******************************'
print*,'Excitation energies (au) (CISD)'
do i = 2, N_states
print*, i ,CI_energy(i) - CI_energy(1)
enddo
print*,''
print*,'******************************'
print*,'Excitation energies (eV) (CISD)'
do i = 2, N_states
print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev
enddo
endif
endif endif
end end

View File

@ -0,0 +1,313 @@
BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_states)]
implicit none
BEGIN_DOC
! One body transition density matrix for all pairs of states n and m, < Psi^n | a_i^\dagger a_a | Psi^m >
END_DOC
integer :: j,k,l,m,k_a,k_b,n
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int)
integer :: exc(0:2,2),n_occ(2)
double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
integer :: krow, kcol, lrow, lcol
PROVIDE psi_det
one_e_tr_dm_mo = 0d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num, &
!$OMP elec_beta_num,one_e_tr_dm_mo,N_det,&
!$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
!$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
!$OMP N_det_alpha_unique,N_det_beta_unique,irp_here)
allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )
tmp_a = 0.d0
!$OMP DO SCHEDULE(dynamic,64)
do k_a=1,N_det
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
! Diagonal part
! -------------
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
do m=1,N_states
do n = 1, N_states
ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,n)
do l=1,elec_alpha_num
j = occ(l,1)
tmp_a(j,j,m,n) += ck
enddo
enddo
enddo
if (k_a == N_det) cycle
l = k_a+1
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
! Fix beta determinant, loop over alphas
do while ( lcol == kcol )
tmp_det2(:) = psi_det_alpha_unique(:, lrow)
call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
if (degree == 1) then
exc = 0
call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
call decode_exc_spin(exc,h1,p1,h2,p2)
do m=1,N_states
do n = 1, N_states
ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,n) * phase
tmp_a(h1,p1,m,n) += ckl
ckl = psi_bilinear_matrix_values(k_a,n)*psi_bilinear_matrix_values(l,m) * phase
tmp_a(p1,h1,m,n) += ckl
enddo
enddo
endif
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_e_tr_dm_mo(:,:,:,:) = one_e_tr_dm_mo(:,:,:,:) + tmp_a(:,:,:,:)
!$OMP END CRITICAL
deallocate(tmp_a)
!$OMP BARRIER
tmp_b = 0.d0
!$OMP DO SCHEDULE(dynamic,64)
do k_b=1,N_det
krow = psi_bilinear_matrix_transp_rows(k_b)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_transp_columns(k_b)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
! Diagonal part
! -------------
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
do m=1,N_states
do n = 1, N_states
ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,n)
do l=1,elec_beta_num
j = occ(l,2)
tmp_b(j,j,m,n) += ck
enddo
enddo
enddo
if (k_b == N_det) cycle
l = k_b+1
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
! Fix beta determinant, loop over alphas
do while ( lrow == krow )
tmp_det2(:) = psi_det_beta_unique(:, lcol)
call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int)
if (degree == 1) then
exc = 0
call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int)
call decode_exc_spin(exc,h1,p1,h2,p2)
do m=1,N_states
do n = 1, N_states
ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,n) * phase
tmp_b(h1,p1,m,n) += ckl
ckl = psi_bilinear_matrix_transp_values(k_b,n)*psi_bilinear_matrix_transp_values(l,m) * phase
tmp_b(p1,h1,m,n) += ckl
enddo
enddo
endif
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_e_tr_dm_mo(:,:,:,:) = one_e_tr_dm_mo(:,:,:,:) + tmp_b(:,:,:,:)
!$OMP END CRITICAL
deallocate(tmp_b)
!$OMP END PARALLEL
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_e_tr_dm_mo_alpha, (mo_num,mo_num,N_states,N_states) ]
&BEGIN_PROVIDER [ double precision, one_e_tr_dm_mo_beta, (mo_num,mo_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! $\alpha$ and $\beta$ one-body transition density matrices for all pairs of states
END_DOC
integer :: j,k,l,m,n,k_a,k_b
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int)
integer :: exc(0:2,2),n_occ(2)
double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
integer :: krow, kcol, lrow, lcol
PROVIDE psi_det
one_e_tr_dm_mo_alpha = 0.d0
one_e_tr_dm_mo_beta = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,k_a,k_b,l,m,n,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num, &
!$OMP elec_beta_num,one_e_tr_dm_mo_alpha,one_e_tr_dm_mo_beta,N_det,&
!$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
!$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
!$OMP N_det_alpha_unique,N_det_beta_unique,irp_here)
allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )
tmp_a = 0.d0
!$OMP DO SCHEDULE(dynamic,64)
do k_a=1,N_det
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
! Diagonal part
! -------------
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
do m=1,N_states
do n = 1, N_states
ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,n)
do l=1,elec_alpha_num
j = occ(l,1)
tmp_a(j,j,m,n) += ck
enddo
enddo
enddo
if (k_a == N_det) cycle
l = k_a+1
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
! Fix beta determinant, loop over alphas
do while ( lcol == kcol )
tmp_det2(:) = psi_det_alpha_unique(:, lrow)
call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
if (degree == 1) then
exc = 0
call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
call decode_exc_spin(exc,h1,p1,h2,p2)
do m=1,N_states
do n = 1, N_states
ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,n) * phase
tmp_a(h1,p1,m,n) += ckl
tmp_a(p1,h1,m,n) += ckl
enddo
enddo
endif
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_e_tr_dm_mo_alpha(:,:,:,:) = one_e_tr_dm_mo_alpha(:,:,:,:) + tmp_a(:,:,:,:)
!$OMP END CRITICAL
deallocate(tmp_a)
tmp_b = 0.d0
!$OMP DO SCHEDULE(dynamic,64)
do k_b=1,N_det
krow = psi_bilinear_matrix_transp_rows(k_b)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_transp_columns(k_b)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
! Diagonal part
! -------------
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
do m=1,N_states
do n = 1, N_states
ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,n)
do l=1,elec_beta_num
j = occ(l,2)
tmp_b(j,j,m,n) += ck
enddo
enddo
enddo
if (k_b == N_det) cycle
l = k_b+1
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
! Fix beta determinant, loop over alphas
do while ( lrow == krow )
tmp_det2(:) = psi_det_beta_unique(:, lcol)
call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int)
if (degree == 1) then
exc = 0
call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int)
call decode_exc_spin(exc,h1,p1,h2,p2)
do m=1,N_states
do n = 1, N_states
ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,n) * phase
tmp_b(h1,p1,m,n) += ckl
tmp_b(p1,h1,m,n) += ckl
enddo
enddo
endif
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_e_tr_dm_mo_beta(:,:,:,:) = one_e_tr_dm_mo_beta(:,:,:,:) + tmp_b(:,:,:,:)
!$OMP END CRITICAL
deallocate(tmp_b)
!$OMP END PARALLEL
END_PROVIDER

View File

@ -0,0 +1,23 @@
[print_all_transitions]
type: logical
doc: If true, print the transition between all the states
interface: ezfio,provider,ocaml
default: false
[calc_dipole_moment]
type: logical
doc: If true, the electric dipole moment will be computed
interface: ezfio,provider,ocaml
default: false
[calc_tr_dipole_moment]
type: logical
doc: If true and N_states > 1, the transition electric dipole moment will be computed
interface: ezfio,provider,ocaml
default: false
[calc_osc_str]
type: logical
doc: If true and N_states > 1, the oscillator strength will be computed
interface: ezfio,provider,ocaml
default: false

2
src/mol_properties/NEED Normal file
View File

@ -0,0 +1,2 @@
determinants
davidson_undressed

View File

@ -0,0 +1,25 @@
# Molecular properties
Available quantities:
- Electric dipole moment
- Electric transition dipole moment
- Oscillator strength
They are not computed by default. To compute them:
```
qp set mol_properties calc_dipole_moment true
qp set mol_properties calc_tr_dipole_moment true
qp set mol_properties calc_osc_str true
```
If you are interested in transitions between two excited states:
```
qp set mol_properties print_all_transitions true
```
They can be obtained by running
```
qp run properties
```
or at each step of a cipsi calculation with
```
qp run fci
```

View File

@ -0,0 +1,13 @@
BEGIN_PROVIDER [double precision, ci_energy_no_diag, (N_states) ]
implicit none
BEGIN_DOC
! CI energy from density matrices and integrals
! Avoid the rediagonalization for ci_energy
END_DOC
ci_energy_no_diag = psi_energy + nuclear_repulsion
END_PROVIDER

View File

@ -0,0 +1,30 @@
BEGIN_PROVIDER [double precision, mo_deriv_1_x , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_deriv_1_y , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_deriv_1_z , (mo_num,mo_num)]
BEGIN_DOC
! array of the integrals of MO_i * d/dx MO_j
! array of the integrals of MO_i * d/dy MO_j
! array of the integrals of MO_i * d/dz MO_j
END_DOC
implicit none
call ao_to_mo( &
ao_deriv_1_x, &
size(ao_deriv_1_x,1), &
mo_deriv_1_x, &
size(mo_deriv_1_x,1) &
)
call ao_to_mo( &
ao_deriv_1_y, &
size(ao_deriv_1_y,1), &
mo_deriv_1_y, &
size(mo_deriv_1_y,1) &
)
call ao_to_mo( &
ao_deriv_1_z, &
size(ao_deriv_1_z,1), &
mo_deriv_1_z, &
size(mo_deriv_1_z,1) &
)
END_PROVIDER

View File

@ -0,0 +1,69 @@
BEGIN_PROVIDER [double precision, multi_s_deriv_1, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_x_deriv_1, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_deriv_1, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_deriv_1, (N_states, N_states)]
implicit none
BEGIN_DOC
! Providers for :
! <Psi_m|v_x|Psi_n>
! <Psi_m|v_y|Psi_n>
! <Psi_m|v_z|Psi_n>
! ||v|| = sqrt(v_x^2 + v_y^2 + v_z^2)
! v_x = d/dx
! Cf. multi_s_dipole_moment for the equations
END_DOC
integer :: istate,jstate ! States
integer :: i,j ! general spatial MOs
double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z
multi_s_x_deriv_1 = 0.d0
multi_s_y_deriv_1 = 0.d0
multi_s_z_deriv_1 = 0.d0
do jstate = 1, N_states
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
multi_s_x_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_x(j,i)
multi_s_y_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_y(j,i)
multi_s_z_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_z(j,i)
enddo
enddo
enddo
enddo
! Nuclei part
nuclei_part_x = 0.d0
nuclei_part_y = 0.d0
nuclei_part_z = 0.d0
do i = 1,nucl_num
nuclei_part_x += nucl_charge(i) * nucl_coord(i,1)
nuclei_part_y += nucl_charge(i) * nucl_coord(i,2)
nuclei_part_z += nucl_charge(i) * nucl_coord(i,3)
enddo
! Only if istate = jstate, otherwise 0 by the orthogonality of the states
do istate = 1, N_states
multi_s_x_deriv_1(istate,istate) += nuclei_part_x
multi_s_y_deriv_1(istate,istate) += nuclei_part_y
multi_s_z_deriv_1(istate,istate) += nuclei_part_z
enddo
! d = <Psi|r|Psi>
do jstate = 1, N_states
do istate = 1, N_states
multi_s_deriv_1(istate,jstate) = &
dsqrt(multi_s_x_deriv_1(istate,jstate)**2 &
+ multi_s_y_deriv_1(istate,jstate)**2 &
+ multi_s_z_deriv_1(istate,jstate)**2)
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,93 @@
! Providers for the dipole moments along x,y,z and the total dipole
! moments.
! The dipole moment along the x axis is:
! \begin{align*}
! \mu_x = < \Psi_m | \sum_i x_i + \sum_A Z_A R_A | \Psi_n >
! \end{align*}
! where $i$ is used for the electrons and $A$ for the nuclei.
! $Z_A$ the charge of the nucleus $A$ and $R_A$ its position in the
! space.
! And it can be computed using the (transition, if n /= m) density
! matrix as a expectation value
! \begin{align*}
! <\Psi_n|x| \Psi_m > = \sum_p \gamma_{pp}^{nm} < \phi_p | x | \phi_p >
! + \sum_{pq, p \neq q} \gamma_{pq}^{nm} < \phi_p | x | \phi_q > + < \Psi_m | \sum_A Z_A R_A | \Psi_n >
! \end{align*}
BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment, (N_states, N_states)]
implicit none
BEGIN_DOC
! Providers for :
! <\Psi_m|\mu_x|\Psi_n>
! <\Psi_m|\mu_y|\Psi_n>
! <\Psi_m|\mu_z|\Psi_n>
! ||\mu|| = \sqrt{\mu_x^2 + \mu_y^2 + \mu_z^2}
!
! <\Psi_n|x| \Psi_m > = \sum_p \gamma_{pp}^{nm} \bra{\phi_p} x \ket{\phi_p}
! + \sum_{pq, p \neq q} \gamma_{pq}^{nm} \bra{\phi_p} x \ket{\phi_q}
! \Psi: wf
! n,m indexes for the states
! p,q: general spatial MOs
! gamma^{nm}: density matrix \bra{\Psi^n} a^{\dagger}_a a_i \ket{\Psi^m}
END_DOC
integer :: istate,jstate ! States
integer :: i,j ! general spatial MOs
double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z
multi_s_x_dipole_moment = 0.d0
multi_s_y_dipole_moment = 0.d0
multi_s_z_dipole_moment = 0.d0
do jstate = 1, N_states
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_x(j,i)
multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_y(j,i)
multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_z(j,i)
enddo
enddo
enddo
enddo
! Nuclei part
nuclei_part_x = 0.d0
nuclei_part_y = 0.d0
nuclei_part_z = 0.d0
do i = 1,nucl_num
nuclei_part_x += nucl_charge(i) * nucl_coord(i,1)
nuclei_part_y += nucl_charge(i) * nucl_coord(i,2)
nuclei_part_z += nucl_charge(i) * nucl_coord(i,3)
enddo
! Only if istate = jstate, otherwise 0 by the orthogonality of the states
do istate = 1, N_states
multi_s_x_dipole_moment(istate,istate) += nuclei_part_x
multi_s_y_dipole_moment(istate,istate) += nuclei_part_y
multi_s_z_dipole_moment(istate,istate) += nuclei_part_z
enddo
! d = <Psi|r|Psi>
do jstate = 1, N_states
do istate = 1, N_states
multi_s_dipole_moment(istate,jstate) = &
dsqrt(multi_s_x_dipole_moment(istate,jstate)**2 &
+ multi_s_y_dipole_moment(istate,jstate)**2 &
+ multi_s_z_dipole_moment(istate,jstate)**2)
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,24 @@
subroutine print_mol_properties()
implicit none
BEGIN_DOC
! Run the propertie calculations
END_DOC
! Electric dipole moment
if (calc_dipole_moment) then
call print_dipole_moment
endif
! Transition electric dipole moment
if (calc_tr_dipole_moment .and. N_states > 1) then
call print_transition_dipole_moment
endif
! Oscillator strength
if (calc_osc_str .and. N_states > 1) then
call print_oscillator_strength
endif
end

View File

@ -0,0 +1,194 @@
! Dipole moments
! Provided
! | N_states | integer | Number of states |
! | multi_s_x_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along x axis |
! | multi_s_y_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along y axis |
! | multi_s_z_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along z axis |
! | multi_s_dipole_moment(N_states,N_states) | double precision | Total (transition) dipole moments |
subroutine print_dipole_moment
implicit none
BEGIN_DOC
! To print the dipole moment ||<\Psi_i|µ|\Psi_i>|| and its x,y,z components
END_DOC
integer :: istate
double precision, allocatable :: d(:), d_x(:), d_y(:), d_z(:)
allocate(d(N_states),d_x(N_states),d_y(N_states),d_z(N_states))
do istate = 1, N_states
d_x(istate) = multi_s_x_dipole_moment(istate,istate)
d_y(istate) = multi_s_y_dipole_moment(istate,istate)
d_z(istate) = multi_s_z_dipole_moment(istate,istate)
d(istate) = multi_s_dipole_moment(istate,istate)
enddo
! Atomic units
print*,''
print*,'# Dipoles:'
print*,'=============================================='
print*,' Dipole moments (au)'
print*,' State X Y Z ||µ||'
do istate = 1, N_states
write(*,'(I5,4(F12.6))') (istate-1), d_x(istate), d_y(istate), d_z(istate), d(istate)
enddo
! Debye
print*,''
print*,' Dipole moments (D)'
print*,' State X Y Z ||µ||'
do istate = 1, N_states
write(*,'(I5,4(F12.6))') (istate-1), d_x(istate)*au_to_D, d_y(istate)*au_to_D, d_z(istate)*au_to_D, d(istate)*au_to_D
enddo
print*,'=============================================='
print*,''
deallocate(d,d_x,d_y,d_z)
end
! Transition dipole moments
! Provided
! | N_states | integer | Number of states |
! | multi_s_x_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along x axis |
! | multi_s_y_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along y axis |
! | multi_s_z_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along z axis |
! | multi_s_dipole_moment(N_states,N_states) | double precision | Total (transition) dipole moments |
subroutine print_transition_dipole_moment
implicit none
BEGIN_DOC
! To print the transition dipole moment ||<\Psi_i|µ|\Psi_j>|| and its components along x, y and z
END_DOC
integer :: istate,jstate, n_states_print
double precision :: f, d, d_x, d_y, d_z, dip_str
if (N_states == 1 .or. N_det == 1) then
return
endif
print*,''
print*,'# Transition dipoles:'
print*,'=============================================='
print*,' Transition dipole moments (au)'
write(*,'(A89)') ' # Transition X Y Z ||µ|| Dip. str. Osc. str.'
if (print_all_transitions) then
n_states_print = N_states
else
n_states_print = 1
endif
do jstate = 1, n_states_print !N_states
do istate = jstate + 1, N_states
d_x = multi_s_x_dipole_moment(istate,jstate)
d_y = multi_s_y_dipole_moment(istate,jstate)
d_z = multi_s_z_dipole_moment(istate,jstate)
dip_str = d_x**2 + d_y**2 + d_z**2
d = multi_s_dipole_moment(istate,jstate)
f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
enddo
enddo
print*,''
print*,' Transition dipole moments (D)'
write(*,'(A89)') ' # Transition X Y Z ||µ|| Dip. str. Osc. str.'
do jstate = 1, n_states_print !N_states
do istate = jstate + 1, N_states
d_x = multi_s_x_dipole_moment(istate,jstate) * au_to_D
d_y = multi_s_y_dipole_moment(istate,jstate) * au_to_D
d_z = multi_s_z_dipole_moment(istate,jstate) * au_to_D
d = multi_s_dipole_moment(istate,jstate)
dip_str = d_x**2 + d_y**2 + d_z**2
f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
d = multi_s_dipole_moment(istate,jstate) * au_to_D
write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
enddo
enddo
print*,'=============================================='
print*,''
end
! Oscillator strengths
! Provided
! | N_states | integer | Number of states |
! | multi_s_dipole_moment(N_states,N_states) | double precision | Total (transition) dipole moments |
! | multi_s_deriv1_moment(N_states,N_states) | double precision | Total (transition) ... |
! | ci_energy_no_diag(N_states) | double precision | CI energy of each state |
! Internal
! | f_l | double precision | Oscillator strength in length gauge |
! | f_v | double precision | Oscillator strength in velocity gauge |
! | f_m | double precision | Oscillator strength in mixed gauge |
! | n_states_print | integer | Number of printed states |
subroutine print_oscillator_strength
implicit none
BEGIN_DOC
! https://doi.org/10.1016/j.cplett.2004.03.126
! Oscillator strength in:
! - length gauge, f^l_{ij} = 2/3 (E_i - E_j) <\Psi_i|r|\Psi_j> <\Psi_j|r|\Psi_i>
! - velocity gauge, f^v_{ij} = 2/3 (E_i - E_j)^(-1) <\Psi_i|v|\Psi_j> <\Psi_j|v|\Psi_i>
! - mixed gauge, f^m_{ij} = -2i/3 <\Psi_i|r|\Psi_j> <\Psi_j|v|\Psi_i>
END_DOC
integer :: istate,jstate,k, n_states_print
double precision :: f_l,f_v,f_m,d,v
if (N_states == 1 .or. N_det == 1) then
return
endif
print*,''
print*,'# Oscillator strength:'
print*,'=============================================='
if (print_all_transitions) then
n_states_print = N_states
else
n_states_print = 1
endif
write(*,'(A103)') ' Oscillator strength in length gauge (f_l), velocity gauge (f_v) and mixed length-velocity gauge (f_m)'
do jstate = 1, n_states_print !N_states
do istate = jstate + 1, N_states
d = multi_s_dipole_moment(istate,jstate)
v = multi_s_deriv_1(istate,jstate)
! Length gauge
f_l = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
! Velocity gauge
f_v = 2d0/3d0 * v * v * 1d0/dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
! Mixed gauge
f_m = 2d0/3d0 * d * v
write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') ' # Transition n.', (istate-1), ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, &
' eV ( ',dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*Ha_to_nm,' nm), f_l=',f_l, ', f_v=', f_v, ', f_m=', f_m, ', <S^2>=', s2_values(istate)
!write(*,'(I4,I4,A4,I3,A6,F6.1,A6,F6.1)') (istate-1), (jstate-1), ' ->', (istate-1), ', %T1=', percent_exc(2,istate), ', %T2=',percent_exc(3,istate)
enddo
enddo
print*,'=============================================='
print*,''
end

View File

@ -0,0 +1,14 @@
program mol_properties
implicit none
BEGIN_DOC
! Calculation of the properties
END_DOC
read_wf = .True.
touch read_wf
call print_mol_properties()
end

View File

@ -1,22 +1,32 @@
BEGIN_PROVIDER [double precision, ha_to_ev] BEGIN_PROVIDER [double precision, ha_to_ev]
&BEGIN_PROVIDER [double precision, au_to_D]
&BEGIN_PROVIDER [double precision, planck_cte]
&BEGIN_PROVIDER [double precision, light_speed]
&BEGIN_PROVIDER [double precision, Ha_to_J]
&BEGIN_PROVIDER [double precision, Ha_to_nm]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Converstion from Hartree to eV ! Some conversion between different units
END_DOC END_DOC
ha_to_ev = 27.211396641308d0 ! Hartree to eV
Ha_to_eV = 27.211396641308d0
END_PROVIDER
BEGIN_PROVIDER [double precision, au_to_D]
implicit none
BEGIN_DOC
! Converstion from au to Debye
END_DOC
! au to Debye
au_to_D = 2.5415802529d0 au_to_D = 2.5415802529d0
END_PROVIDER ! Planck's constant in SI units
planck_cte = 6.62606957d-34
! Light speed in SI units
light_speed = 2.99792458d10
! Hartree to Joule
Ha_to_J = 4.35974434d-18
! Hartree to nm
Ha_to_nm = 1d9 * (planck_cte * light_speed) / Ha_to_J
END_PROVIDER

77
src/utils_cc/EZFIO.cfg Normal file
View File

@ -0,0 +1,77 @@
[cc_thresh_conv]
type: double precision
doc: Threshold for the convergence of the residual equations.
interface: ezfio,ocaml,provider
default: 1e-6
[cc_max_iter]
type: integer
doc: Maximum number of iterations.
interface: ezfio,ocaml,provider
default: 100
[cc_diis_depth]
type: integer
doc: Maximum depth of the DIIS, i.e., maximum number of iterations that the DIIS keeps in memory. Warning, we allocate matrices with the diis depth at the beginning without update. If you don't have enough memory it should crash in memory.
interface: ezfio,ocaml,provider
default: 8
[cc_level_shift]
type: double precision
doc: Level shift for the CC
interface: ezfio,ocaml,provider
default: 0.0
[cc_level_shift_guess]
type: double precision
doc: Level shift for the guess of the CC amplitudes
interface: ezfio,ocaml,provider
default: 0.0
[cc_update_method]
type: character*(32)
doc: Method used to update the CC amplitudes. none -> normal, diis -> with diis.
interface: ezfio,ocaml,provider
default: diis
[cc_guess_t1]
type: character*(32)
doc: Guess used to initialize the T1 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk.
interface: ezfio,ocaml,provider
default: MP
[cc_guess_t2]
type: character*(32)
doc: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk.
interface: ezfio,ocaml,provider
default: MP
[cc_write_t1]
type: logical
doc: If true, it will write on disk the T1 amplitudes at the end of the calculation.
interface: ezfio,ocaml,provider
default: False
[cc_write_t2]
type: logical
doc: If true, it will write on disk the T2 amplitudes at the end of the calculation.
interface: ezfio,ocaml,provider
default: False
[cc_par_t]
type: logical
doc: If true, the CCSD(T) will be computed.
interface: ezfio,ocaml,provider
default: False
[cc_dev]
type: logical
doc: Only for dev purposes.
interface: ezfio,ocaml,provider
default: False
[cc_ref]
type: integer
doc: Index of the reference determinant in psi_det for CC calculation.
interface: ezfio,ocaml,provider
default: 1

4
src/utils_cc/NEED Normal file
View File

@ -0,0 +1,4 @@
hartree_fock
two_body_rdm
bitmask
determinants

34
src/utils_cc/README.md Normal file
View File

@ -0,0 +1,34 @@
# Utils for CC
Utils for the CC modules.
## Contents
- Providers related to reference occupancy
- Integrals related to the reference
- Diis for CC (but can be used for something else if you provide your own error vector)
- Guess for CC amplitudes
- Routines to update the CC amplitudes
- Phase between to arbitrary determinants
- print of the qp edit wf
## Keywords
- cc_thresh_conv: Threshold for the convergence of the residual equations. Default: 1e-6.
- cc_max_iter: Maximum number of iterations. Default: 100.
- cc_diis_depth: Diis depth. Default: 8.
- cc_level_shift: Level shift for the CC. Default: 0.0.
- cc_level_shift_guess: Level shift for the MP guess of the amplitudes. Default: 0.0.
- cc_update_method: Method used to update the CC amplitudes. none -> normal, diis -> with diis. Default: diis.
- cc_guess_t1: Guess used to initialize the T1 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk. Default: MP.
- cc_guess_t2: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk. Default: MP.
- cc_write_t1: If true, it will write on disk the T1 amplitudes at the end of the calculation. Default: False.
- cc_write_t2: If true, it will write on disk the T2 amplitudes at the end of the calculation. Default: False.
- cc_par_t: If true, the CCSD(T) will be computed.
- cc_ref: Index of the reference determinant in psi_det for CC calculation. Default: 1.
## Org files
The org files are stored in the directory org in order to avoid overwriting on user changes.
The org files can be modified, to export the change to the source code, run
```
./TANGLE_org_mode.sh and
mv *.irp.f ../.
```

529
src/utils_cc/diis.irp.f Normal file
View File

@ -0,0 +1,529 @@
! Code
subroutine diis_cc(all_err,all_t,sze,m,iter,t)
implicit none
BEGIN_DOC
! DIIS. Take the error vectors and the amplitudes of the previous
! iterations to compute the new amplitudes
END_DOC
! {err_i}_{i=1}^{m_it} -> B -> c
! {t_i}_{i=1}^{m_it}, c, {err_i}_{i=1}^{m_it} -> t_{m_it+1}
integer, intent(in) :: m,iter,sze
double precision, intent(in) :: all_err(sze,m)
double precision, intent(in) :: all_t(sze,m)
double precision, intent(out) :: t(sze)
double precision, allocatable :: B(:,:), c(:), zero(:)
integer :: m_iter
integer :: i,j,k
integer :: info
integer, allocatable :: ipiv(:)
double precision :: accu
m_iter = min(m,iter)
!print*,'m_iter',m_iter
allocate(B(m_iter+1,m_iter+1), c(m_iter), zero(m_iter+1))
allocate(ipiv(m+1))
! B(i,j) = < err(iter-m_iter+j),err(iter-m_iter+i) > ! iter-m_iter will be zero for us
B = 0d0
!$OMP PARALLEL &
!$OMP SHARED(B,m,m_iter,sze,all_err) &
!$OMP PRIVATE(i,j,k,accu) &
!$OMP DEFAULT(NONE)
do j = 1, m_iter
do i = 1, m_iter
accu = 0d0
!$OMP DO
do k = 1, sze
! the errors of the ith iteration are in all_err(:,m+1-i)
accu = accu + all_err(k,m+1-i) * all_err(k,m+1-j)
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
B(i,j) = B(i,j) + accu
!$OMP END CRITICAL
enddo
enddo
!$OMP END PARALLEL
do i = 1, m_iter
B(i,m_iter+1) = -1
enddo
do j = 1, m_iter
B(m_iter+1,j) = -1
enddo
! Debug
!print*,'B'
!do i = 1, m_iter+1
! write(*,'(100(F10.6))') B(i,:)
!enddo
! (0 0 .... 0 -1)
zero = 0d0
zero(m_iter+1) = -1d0
! Solve B.c = zero
call dgesv(m_iter+1, 1, B, size(B,1), ipiv, zero, size(zero,1), info)
if (info /= 0) then
print*,'DIIS error in dgesv:', info
call abort
endif
! c corresponds to the m_iter first solutions
c = zero(1:m_iter)
! Debug
!print*,'c',c
!print*,'all_t'
!do i = 1, m
! write(*,'(100(F10.6))') all_t(:,i)
!enddo
!print*,'all_err'
!do i = 1, m
! write(*,'(100(F10.6))') all_err(:,i)
!enddo
! update T
!$OMP PARALLEL &
!$OMP SHARED(t,c,m,all_err,all_t,sze,m_iter) &
!$OMP PRIVATE(i,j,accu) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, sze
t(i) = 0d0
enddo
!$OMP END DO
do i = 1, m_iter
!$OMP DO
do j = 1, sze
t(j) = t(j) + c(i) * (all_t(j,m+1-i) + all_err(j,m+1-i))
enddo
!$OMP END DO
enddo
!$OMP END PARALLEL
!print*,'new t',t
deallocate(ipiv,B,c,zero)
end
! Update all err
subroutine update_all_err(err,all_err,sze,m,iter)
implicit none
BEGIN_DOC
! Shift all the err vectors of the previous iterations to add the new one
! The last err vector is placed in the last position and all the others are
! moved toward the first one.
END_DOC
integer, intent(in) :: m, iter, sze
double precision, intent(in) :: err(sze)
double precision, intent(inout) :: all_err(sze,m)
integer :: i,j
integer :: m_iter
m_iter = min(m,iter)
! Shift
!$OMP PARALLEL &
!$OMP SHARED(m,all_err,err,sze) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do i = 1, m-1
!$OMP DO
do j = 1, sze
all_err(j,i) = all_err(j,i+1)
enddo
!$OMP END DO
enddo
! Debug
!print*,'shift err'
!do i = 1, m
! print*,i, all_err(:,i)
!enddo
! New
!$OMP DO
do i = 1, sze
all_err(i,m) = err(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Debug
!print*,'Updated err'
!do i = 1, m
! print*,i, all_err(:,i)
!enddo
end
! Update all t
subroutine update_all_t(t,all_t,sze,m,iter)
implicit none
BEGIN_DOC
! Shift all the t vectors of the previous iterations to add the new one
! The last t vector is placed in the last position and all the others are
! moved toward the first one.
END_DOC
integer, intent(in) :: m, iter, sze
double precision, intent(in) :: t(sze)
double precision, intent(inout) :: all_t(sze,m)
integer :: i,j
integer :: m_iter
m_iter = min(m,iter)
! Shift
!$OMP PARALLEL &
!$OMP SHARED(m,all_t,t,sze) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do i = 1, m-1
!$OMP DO
do j = 1, sze
all_t(j,i) = all_t(j,i+1)
enddo
!$OMP END DO
enddo
! New
!$OMP DO
do i = 1, sze
all_t(i,m) = t(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Debug
!print*,'Updated t'
!do i = 1, m
! print*,i, all_t(:,i)
!enddo
end
! Err1
subroutine compute_err1(nO,nV,f_o,f_v,r1,err1)
implicit none
BEGIN_DOC
! Compute the error vector for the t1
END_DOC
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r1(nO,nV)
double precision, intent(out) :: err1(nO,nV)
integer :: i,a
!$OMP PARALLEL &
!$OMP SHARED(err1,r1,f_o,f_v,nO,nV,cc_level_shift) &
!$OMP PRIVATE(i,a) &
!$OMP DEFAULT(NONE)
!$OMP DO
do a = 1, nV
do i = 1, nO
err1(i,a) = - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
! Err2
subroutine compute_err2(nO,nV,f_o,f_v,r2,err2)
implicit none
BEGIN_DOC
! Compute the error vector for the t2
END_DOC
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r2(nO,nO,nV,nV)
double precision, intent(out) :: err2(nO,nO,nV,nV)
integer :: i,j,a,b
!$OMP PARALLEL &
!$OMP SHARED(err2,r2,f_o,f_v,nO,nV,cc_level_shift) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
err2(i,j,a,b) = - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
! Update t
subroutine update_t_ccsd(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: err1(:,:), err2(:,:,:,:)
double precision, allocatable :: tmp_err1(:), tmp_err2(:)
double precision, allocatable :: tmp_t1(:), tmp_t2(:)
if (cc_update_method == 'diis') then
allocate(err1(nO,nV), err2(nO,nO,nV,nV))
allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
! DIIS T1, it is not always good since the t1 can be small
! That's why there is a call to update the t1 in the standard way
! T1 error tensor
!call compute_err1(nO,nV,f_o,f_v,r1,err1)
! Transfo errors and parameters in vectors
!tmp_err1 = reshape(err1,(/nO*nV/))
!tmp_t1 = reshape(t1 ,(/nO*nV/))
! Add the error and parameter vectors with those of the previous iterations
!call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
!call update_all_t (tmp_t1 ,all_t1 ,nO*nV,cc_diis_depth,nb_iter+1)
! Diis and reshape T as a tensor
!call diis_cc(all_err1,all_t1,nO*nV,cc_diis_depth,nb_iter+1,tmp_t1)
!t1 = reshape(tmp_t1 ,(/nO,nV/))
call update_t1(nO,nV,f_o,f_v,r1,t1)
! DIIS T2
! T2 error tensor
call compute_err2(nO,nV,f_o,f_v,r2,err2)
! Transfo errors and parameters in vectors
tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
tmp_t2 = reshape(t2 ,(/nO*nO*nV*nV/))
! Add the error and parameter vectors with those of the previous iterations
call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t2 ,all_t2 ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Diis and reshape T as a tensor
call diis_cc(all_err2,all_t2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t2)
t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2)
! Standard update as T = T - Delta
elseif (cc_update_method == 'none') then
call update_t1(nO,nV,f_o,f_v,r1,t1)
call update_t2(nO,nV,f_o,f_v,r2,t2)
else
print*,'Unkonw cc_method_method: '//cc_update_method
endif
end
! Update t v2
subroutine update_t_ccsd_diis(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: all_t(:,:), all_err(:,:), tmp_t(:)
double precision, allocatable :: err1(:,:), err2(:,:,:,:)
double precision, allocatable :: tmp_err1(:), tmp_err2(:)
double precision, allocatable :: tmp_t1(:), tmp_t2(:)
integer :: i,j
! Allocate
allocate(all_err(nO*nV+nO*nO*nV*nV,cc_diis_depth), all_t(nO*nV+nO*nO*nV*nV,cc_diis_depth))
allocate(tmp_t(nO*nV+nO*nO*nV*nV))
allocate(err1(nO,nV), err2(nO,nO,nV,nV))
allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
! Compute the errors and reshape them as vector
call compute_err1(nO,nV,f_o,f_v,r1,err1)
call compute_err2(nO,nV,f_o,f_v,r2,err2)
tmp_err1 = reshape(err1,(/nO*nV/))
tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
tmp_t1 = reshape(t1 ,(/nO*nV/))
tmp_t2 = reshape(t2 ,(/nO*nO*nV*nV/))
! Update the errors and parameters for the diis
call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t1 ,all_t1 ,nO*nV,cc_diis_depth,nb_iter+1)
call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t2 ,all_t2 ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Gather the different parameters and errors
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,all_err,all_err1,all_err2,cc_diis_depth,&
!$OMP all_t,all_t1,all_t2) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nV
all_err(i,j) = all_err1(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nO*nV*nV
all_err(i+nO*nV,j) = all_err2(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nV
all_t(i,j) = all_t1(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nO*nV*nV
all_t(i+nO*nV,j) = all_t2(i,j)
enddo
!$OMP END DO
enddo
!$OMP END PARALLEL
! Diis
call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t)
! Split the resulting vector
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp_t,tmp_t1,tmp_t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
tmp_t1(i) = tmp_t(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
tmp_t2(i) = tmp_t(i+nO*nV)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Reshape as tensors
t1 = reshape(tmp_t1 ,(/nO,nV/))
t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
! Deallocate
deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2,all_t,all_err)
end
! Update t v3
subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO*nV), t2(nO*nO*nV*nV)
double precision, intent(inout) :: all_err(nO*nV+nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t(nO*nV+nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: tmp(:)
integer :: i,j
! Allocate
allocate(tmp(nO*nV+nO*nO*nV*nV))
! Compute the errors
call compute_err1(nO,nV,f_o,f_v,r1,tmp(1:nO*nV))
call compute_err2(nO,nV,f_o,f_v,r2,tmp(nO*nV+1:nO*nV+nO*nO*nV*nV))
! Update the errors and parameters for the diis
call update_all_err(tmp,all_err,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp,t1,t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
tmp(i) = t1(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
tmp(i+nO*nV) = t2(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
call update_all_t(tmp,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Diis
call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp)
! Split the resulting vector
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp,t1,t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
t1(i) = tmp(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
t2(i) = tmp(i+nO*nV)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Deallocate
deallocate(tmp)
end

13
src/utils_cc/energy.irp.f Normal file
View File

@ -0,0 +1,13 @@
subroutine det_energy(det,energy)
implicit none
integer(bit_kind), intent(in) :: det
double precision, intent(out) :: energy
call i_H_j(det,det,N_int,energy)
energy = energy + nuclear_repulsion
end

213
src/utils_cc/guess_t.irp.f Normal file
View File

@ -0,0 +1,213 @@
! T1
subroutine guess_t1(nO,nV,f_o,f_v,f_ov,t1)
implicit none
BEGIN_DOC
! Update the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
! inout
double precision, intent(out) :: t1(nO, nV)
! internal
integer :: i,a
if (trim(cc_guess_t1) == 'none') then
t1 = 0d0
else if (trim(cc_guess_t1) == 'MP') then
do a = 1, nV
do i = 1, nO
t1(i,a) = f_ov(i,a) / (f_o(i) - f_v(a) - cc_level_shift_guess)
enddo
enddo
else if (trim(cc_guess_t1) == 'read') then
call read_t1(nO,nV,t1)
else
print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t1)
call abort
endif
end
! T2
subroutine guess_t2(nO,nV,f_o,f_v,v_oovv,t2)
implicit none
BEGIN_DOC
! Update the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), v_oovv(nO, nO, nV, nV)
! inout
double precision, intent(out) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
if (trim(cc_guess_t2) == 'none') then
t2 = 0d0
else if (trim(cc_guess_t2) == 'MP') then
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
t2(i,j,a,b) = v_oovv(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift_guess)
enddo
enddo
enddo
enddo
else if (trim(cc_guess_t2) == 'read') then
call read_t2(nO,nV,t2)
else
print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t2)
call abort
endif
end
! T1
subroutine write_t1(nO,nV,t1)
implicit none
BEGIN_DOC
! Write the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO, nV)
! internal
integer :: i,a
if (cc_write_t1) then
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
do a = 1, nV
do i = 1, nO
write(11,'(F20.12)') t1(i,a)
enddo
enddo
close(11)
endif
end
! T2
subroutine write_t2(nO,nV,t2)
implicit none
BEGIN_DOC
! Write the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
if (cc_write_t2) then
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
write(11,'(F20.12)') t2(i,j,a,b)
enddo
enddo
enddo
enddo
close(11)
endif
end
! T1
subroutine read_t1(nO,nV,t1)
implicit none
BEGIN_DOC
! Read the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(out) :: t1(nO, nV)
! internal
integer :: i,a
logical :: ok
inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
if (.not. ok) then
print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
print*, 'Do a first calculation with cc_write_t1 = True'
print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
call abort
endif
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
do a = 1, nV
do i = 1, nO
read(11,'(F20.12)') t1(i,a)
enddo
enddo
close(11)
end
! T2
subroutine read_t2(nO,nV,t2)
implicit none
BEGIN_DOC
! Read the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(out) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
logical :: ok
inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
if (.not. ok) then
print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
print*, 'Do a first calculation with cc_write_t2 = True'
print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
call abort
endif
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
read(11,'(F20.12)') t2(i,j,a,b)
enddo
enddo
enddo
enddo
close(11)
end

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,328 @@
! N spin orb
subroutine extract_n_spin(det,n)
implicit none
BEGIN_DOC
! Returns the number of occupied alpha, occupied beta, virtual alpha, virtual beta spin orbitals
! in det without counting the core and deleted orbitals in the format n(nOa,nOb,nVa,nVb)
END_DOC
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: n(4)
integer(bit_kind) :: res(N_int,2)
integer :: i, si
logical :: ok, is_core, is_del
! Init
n = 0
! Loop over the spin
do si = 1, 2
do i = 1, mo_num
call apply_hole(det, si, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
n(si) = n(si) + 1
else
! hole
n(si+2) = n(si+2) + 1
endif
enddo
enddo
!print*,n(1),n(2),n(3),n(4)
end
! Spin
subroutine extract_list_orb_spin(det,nO_m,nV_m,list_occ,list_vir)
implicit none
BEGIN_DOC
! Returns the the list of occupied alpha/beta, virtual alpha/beta spin orbitals
! size(nO_m,1) must be max(nOa,nOb) and size(nV_m,1) must be max(nVa,nVb)
END_DOC
integer, intent(in) :: nO_m, nV_m
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: list_occ(nO_m,2), list_vir(nV_m,2)
integer(bit_kind) :: res(N_int,2)
integer :: i, si, idx_o, idx_v, idx_i, idx_b
logical :: ok, is_core, is_del
list_occ = 0
list_vir = 0
! List of occ/vir alpha/beta
! occ alpha -> list_occ(:,1)
! occ beta -> list_occ(:,2)
! vir alpha -> list_vir(:,1)
! vir beta -> list_vir(:,2)
! Loop over the spin
do si = 1, 2
! tmp idx
idx_o = 1
idx_v = 1
do i = 1, mo_num
call apply_hole(det, si, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
list_occ(idx_o,si) = i
idx_o = idx_o + 1
else
! hole
list_vir(idx_v,si) = i
idx_v = idx_v + 1
endif
enddo
enddo
end
! Space
subroutine extract_list_orb_space(det,nO,nV,list_occ,list_vir)
implicit none
BEGIN_DOC
! Returns the the list of occupied and virtual alpha spin orbitals
END_DOC
integer, intent(in) :: nO, nV
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: list_occ(nO), list_vir(nV)
integer(bit_kind) :: res(N_int,2)
integer :: i, si, idx_o, idx_v, idx_i, idx_b
logical :: ok, is_core, is_del
if (elec_alpha_num /= elec_beta_num) then
print*,'Error elec_alpha_num /= elec_beta_num, impossible to create cc_list_occ and cc_list_vir, abort'
call abort
endif
list_occ = 0
list_vir = 0
! List of occ/vir alpha
! occ alpha -> list_occ(:,1)
! vir alpha -> list_vir(:,1)
! tmp idx
idx_o = 1
idx_v = 1
do i = 1, mo_num
call apply_hole(det, 1, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
list_occ(idx_o) = i
idx_o = idx_o + 1
else
! hole
list_vir(idx_v) = i
idx_v = idx_v + 1
endif
enddo
end
! is_core
function is_core(i)
implicit none
BEGIN_DOC
! True if the orbital i is a core orbital
END_DOC
integer, intent(in) :: i
logical :: is_core
integer :: j
! Init
is_core = .False.
! Search
do j = 1, dim_list_core_orb
if (list_core(j) == i) then
is_core = .True.
exit
endif
enddo
end
! is_del
function is_del(i)
implicit none
BEGIN_DOC
! True if the orbital i is a deleted orbital
END_DOC
integer, intent(in) :: i
logical :: is_del
integer :: j
! Init
is_del = .False.
! Search
do j = 1, dim_list_core_orb
if (list_core(j) == i) then
is_del = .True.
exit
endif
enddo
end
! N orb
BEGIN_PROVIDER [integer, cc_nO_m]
&BEGIN_PROVIDER [integer, cc_nOa]
&BEGIN_PROVIDER [integer, cc_nOb]
&BEGIN_PROVIDER [integer, cc_nOab]
&BEGIN_PROVIDER [integer, cc_nV_m]
&BEGIN_PROVIDER [integer, cc_nVa]
&BEGIN_PROVIDER [integer, cc_nVb]
&BEGIN_PROVIDER [integer, cc_nVab]
&BEGIN_PROVIDER [integer, cc_n_mo]
&BEGIN_PROVIDER [integer, cc_nO_S, (2)]
&BEGIN_PROVIDER [integer, cc_nV_S, (2)]
implicit none
BEGIN_DOC
! Number of orbitals without core and deleted ones of the cc_ref det in psi_det
! a: alpha, b: beta
! nO_m: max(a,b) occupied
! nOa: nb a occupied
! nOb: nb b occupied
! nOab: nb a+b occupied
! nV_m: max(a,b) virtual
! nVa: nb a virtual
! nVb: nb b virtual
! nVab: nb a+b virtual
END_DOC
integer :: n_spin(4)
! Extract number of occ/vir alpha/beta spin orbitals
call extract_n_spin(psi_det(1,1,cc_ref),n_spin)
cc_nOa = n_spin(1)
cc_nOb = n_spin(2)
cc_nOab = cc_nOa + cc_nOb !n_spin(1) + n_spin(2)
cc_nO_m = max(cc_nOa,cc_nOb) !max(n_spin(1), n_spin(2))
cc_nVa = n_spin(3)
cc_nVb = n_spin(4)
cc_nVab = cc_nVa + cc_nVb !n_spin(3) + n_spin(4)
cc_nV_m = max(cc_nVa,cc_nVb) !max(n_spin(3), n_spin(4))
cc_n_mo = cc_nVa + cc_nVb !n_spin(1) + n_spin(3)
cc_nO_S = (/cc_nOa,cc_nOb/)
cc_nV_S = (/cc_nVa,cc_nVb/)
END_PROVIDER
! General
BEGIN_PROVIDER [integer, cc_list_gen, (cc_n_mo)]
implicit none
BEGIN_DOC
! List of general orbitals without core and deleted ones
END_DOC
integer :: i,j
logical :: is_core, is_del
j = 1
do i = 1, mo_num
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
cc_list_gen(j) = i
j = j+1
enddo
END_PROVIDER
! Space
BEGIN_PROVIDER [integer, cc_list_occ, (cc_nOa)]
&BEGIN_PROVIDER [integer, cc_list_vir, (cc_nVa)]
implicit none
BEGIN_DOC
! List of occupied and virtual spatial orbitals without core and deleted ones
END_DOC
call extract_list_orb_space(psi_det(1,1,cc_ref),cc_nOa,cc_nVa,cc_list_occ,cc_list_vir)
END_PROVIDER
! Spin
BEGIN_PROVIDER [integer, cc_list_occ_spin, (cc_nO_m,2)]
&BEGIN_PROVIDER [integer, cc_list_vir_spin, (cc_nV_m,2)]
&BEGIN_PROVIDER [logical, cc_ref_is_open_shell]
implicit none
BEGIN_DOC
! List of occupied and virtual spin orbitals without core and deleted ones
END_DOC
integer :: i
call extract_list_orb_spin(psi_det(1,1,cc_ref),cc_nO_m,cc_nV_m,cc_list_occ_spin,cc_list_vir_spin)
cc_ref_is_open_shell = .False.
do i = 1, cc_nO_m
if (cc_list_occ_spin(i,1) /= cc_list_occ_spin(i,2)) then
cc_ref_is_open_shell = .True.
endif
enddo
END_PROVIDER

View File

@ -0,0 +1,7 @@
#!/bin/sh
list='ls *.org'
for element in $list
do
emacs --batch $element -f org-babel-tangle
done

574
src/utils_cc/org/diis.org Normal file
View File

@ -0,0 +1,574 @@
* DIIS
https://hal.archives-ouvertes.fr/hal-02492983/document
Maxime Chupin, Mi-Song Dupuy, Guillaume Legendre, Eric Séré. Convergence analysis of adaptive
DIIS algorithms witerh application to electronic ground state calculations.
ESAIM: Mathematical Modelling and Numerical Analysis, EDP Sciences, 2021, 55 (6), pp.2785 - 2825. 10.1051/m2an/2021069ff.ffhal-02492983v5
t_{k+1} = g(t_k)
err_k = f(t_k) = t_{k+1} - t_k
m_k = min(m,k)
m maximal depth
t_{k+1} = \sum_{i=0}^{m_k} c_i^k g(t_{k-m_k+i})
\sum_{i=0}^{m_k} c_i^k = 1
b_{ij}^k = < err^{k-m_k+j}, err^{k-m_k+i} >
(b -1) ( c^k ) = ( 0 )
(-1 0) ( \lambda) ( -1 )
lambda is used to put the constraint \sum_{i=0}^{m_k} c_i^k = 1
In: t_0, err_0, m
err_0 = g(t_0)
k = 0
m_k = 0
while ||err_k|| > CC
A.x=b
t_{k+1} = \sum_{i=0}^{m_k} c_i^k g(t_{k-m_k+i})
err_{k+1} = f(t_{k+1})
m_{k+1} = min(m_k+1,m)
k = k +1
end
* Code
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine diis_cc(all_err,all_t,sze,m,iter,t)
implicit none
BEGIN_DOC
! DIIS. Take the error vectors and the amplitudes of the previous
! iterations to compute the new amplitudes
END_DOC
! {err_i}_{i=1}^{m_it} -> B -> c
! {t_i}_{i=1}^{m_it}, c, {err_i}_{i=1}^{m_it} -> t_{m_it+1}
integer, intent(in) :: m,iter,sze
double precision, intent(in) :: all_err(sze,m)
double precision, intent(in) :: all_t(sze,m)
double precision, intent(out) :: t(sze)
double precision, allocatable :: B(:,:), c(:), zero(:)
integer :: m_iter
integer :: i,j,k
integer :: info
integer, allocatable :: ipiv(:)
double precision :: accu
m_iter = min(m,iter)
!print*,'m_iter',m_iter
allocate(B(m_iter+1,m_iter+1), c(m_iter), zero(m_iter+1))
allocate(ipiv(m+1))
! B(i,j) = < err(iter-m_iter+j),err(iter-m_iter+i) > ! iter-m_iter will be zero for us
B = 0d0
!$OMP PARALLEL &
!$OMP SHARED(B,m,m_iter,sze,all_err) &
!$OMP PRIVATE(i,j,k,accu) &
!$OMP DEFAULT(NONE)
do j = 1, m_iter
do i = 1, m_iter
accu = 0d0
!$OMP DO
do k = 1, sze
! the errors of the ith iteration are in all_err(:,m+1-i)
accu = accu + all_err(k,m+1-i) * all_err(k,m+1-j)
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
B(i,j) = B(i,j) + accu
!$OMP END CRITICAL
enddo
enddo
!$OMP END PARALLEL
do i = 1, m_iter
B(i,m_iter+1) = -1
enddo
do j = 1, m_iter
B(m_iter+1,j) = -1
enddo
! Debug
!print*,'B'
!do i = 1, m_iter+1
! write(*,'(100(F10.6))') B(i,:)
!enddo
! (0 0 .... 0 -1)
zero = 0d0
zero(m_iter+1) = -1d0
! Solve B.c = zero
call dgesv(m_iter+1, 1, B, size(B,1), ipiv, zero, size(zero,1), info)
if (info /= 0) then
print*,'DIIS error in dgesv:', info
call abort
endif
! c corresponds to the m_iter first solutions
c = zero(1:m_iter)
! Debug
!print*,'c',c
!print*,'all_t'
!do i = 1, m
! write(*,'(100(F10.6))') all_t(:,i)
!enddo
!print*,'all_err'
!do i = 1, m
! write(*,'(100(F10.6))') all_err(:,i)
!enddo
! update T
!$OMP PARALLEL &
!$OMP SHARED(t,c,m,all_err,all_t,sze,m_iter) &
!$OMP PRIVATE(i,j,accu) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, sze
t(i) = 0d0
enddo
!$OMP END DO
do i = 1, m_iter
!$OMP DO
do j = 1, sze
t(j) = t(j) + c(i) * (all_t(j,m+1-i) + all_err(j,m+1-i))
enddo
!$OMP END DO
enddo
!$OMP END PARALLEL
!print*,'new t',t
deallocate(ipiv,B,c,zero)
end
#+end_src
** Update all err
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine update_all_err(err,all_err,sze,m,iter)
implicit none
BEGIN_DOC
! Shift all the err vectors of the previous iterations to add the new one
! The last err vector is placed in the last position and all the others are
! moved toward the first one.
END_DOC
integer, intent(in) :: m, iter, sze
double precision, intent(in) :: err(sze)
double precision, intent(inout) :: all_err(sze,m)
integer :: i,j
integer :: m_iter
m_iter = min(m,iter)
! Shift
!$OMP PARALLEL &
!$OMP SHARED(m,all_err,err,sze) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do i = 1, m-1
!$OMP DO
do j = 1, sze
all_err(j,i) = all_err(j,i+1)
enddo
!$OMP END DO
enddo
! Debug
!print*,'shift err'
!do i = 1, m
! print*,i, all_err(:,i)
!enddo
! New
!$OMP DO
do i = 1, sze
all_err(i,m) = err(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Debug
!print*,'Updated err'
!do i = 1, m
! print*,i, all_err(:,i)
!enddo
end
#+end_src
** Update all t
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine update_all_t(t,all_t,sze,m,iter)
implicit none
BEGIN_DOC
! Shift all the t vectors of the previous iterations to add the new one
! The last t vector is placed in the last position and all the others are
! moved toward the first one.
END_DOC
integer, intent(in) :: m, iter, sze
double precision, intent(in) :: t(sze)
double precision, intent(inout) :: all_t(sze,m)
integer :: i,j
integer :: m_iter
m_iter = min(m,iter)
! Shift
!$OMP PARALLEL &
!$OMP SHARED(m,all_t,t,sze) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do i = 1, m-1
!$OMP DO
do j = 1, sze
all_t(j,i) = all_t(j,i+1)
enddo
!$OMP END DO
enddo
! New
!$OMP DO
do i = 1, sze
all_t(i,m) = t(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Debug
!print*,'Updated t'
!do i = 1, m
! print*,i, all_t(:,i)
!enddo
end
#+end_src
** Err
*** Err1
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine compute_err1(nO,nV,f_o,f_v,r1,err1)
implicit none
BEGIN_DOC
! Compute the error vector for the t1
END_DOC
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r1(nO,nV)
double precision, intent(out) :: err1(nO,nV)
integer :: i,a
!$OMP PARALLEL &
!$OMP SHARED(err1,r1,f_o,f_v,nO,nV,cc_level_shift) &
!$OMP PRIVATE(i,a) &
!$OMP DEFAULT(NONE)
!$OMP DO
do a = 1, nV
do i = 1, nO
err1(i,a) = - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
#+end_src
*** Err2
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine compute_err2(nO,nV,f_o,f_v,r2,err2)
implicit none
BEGIN_DOC
! Compute the error vector for the t2
END_DOC
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r2(nO,nO,nV,nV)
double precision, intent(out) :: err2(nO,nO,nV,nV)
integer :: i,j,a,b
!$OMP PARALLEL &
!$OMP SHARED(err2,r2,f_o,f_v,nO,nV,cc_level_shift) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
err2(i,j,a,b) = - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
#+end_src
* Gather call diis
** Update t
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine update_t_ccsd(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: err1(:,:), err2(:,:,:,:)
double precision, allocatable :: tmp_err1(:), tmp_err2(:)
double precision, allocatable :: tmp_t1(:), tmp_t2(:)
if (cc_update_method == 'diis') then
allocate(err1(nO,nV), err2(nO,nO,nV,nV))
allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
! DIIS T1, it is not always good since the t1 can be small
! That's why there is a call to update the t1 in the standard way
! T1 error tensor
!call compute_err1(nO,nV,f_o,f_v,r1,err1)
! Transfo errors and parameters in vectors
!tmp_err1 = reshape(err1,(/nO*nV/))
!tmp_t1 = reshape(t1 ,(/nO*nV/))
! Add the error and parameter vectors with those of the previous iterations
!call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
!call update_all_t (tmp_t1 ,all_t1 ,nO*nV,cc_diis_depth,nb_iter+1)
! Diis and reshape T as a tensor
!call diis_cc(all_err1,all_t1,nO*nV,cc_diis_depth,nb_iter+1,tmp_t1)
!t1 = reshape(tmp_t1 ,(/nO,nV/))
call update_t1(nO,nV,f_o,f_v,r1,t1)
! DIIS T2
! T2 error tensor
call compute_err2(nO,nV,f_o,f_v,r2,err2)
! Transfo errors and parameters in vectors
tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
tmp_t2 = reshape(t2 ,(/nO*nO*nV*nV/))
! Add the error and parameter vectors with those of the previous iterations
call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t2 ,all_t2 ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Diis and reshape T as a tensor
call diis_cc(all_err2,all_t2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t2)
t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2)
! Standard update as T = T - Delta
elseif (cc_update_method == 'none') then
call update_t1(nO,nV,f_o,f_v,r1,t1)
call update_t2(nO,nV,f_o,f_v,r2,t2)
else
print*,'Unkonw cc_method_method: '//cc_update_method
endif
end
#+end_src
** Update t v2
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine update_t_ccsd_diis(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: all_t(:,:), all_err(:,:), tmp_t(:)
double precision, allocatable :: err1(:,:), err2(:,:,:,:)
double precision, allocatable :: tmp_err1(:), tmp_err2(:)
double precision, allocatable :: tmp_t1(:), tmp_t2(:)
integer :: i,j
! Allocate
allocate(all_err(nO*nV+nO*nO*nV*nV,cc_diis_depth), all_t(nO*nV+nO*nO*nV*nV,cc_diis_depth))
allocate(tmp_t(nO*nV+nO*nO*nV*nV))
allocate(err1(nO,nV), err2(nO,nO,nV,nV))
allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
! Compute the errors and reshape them as vector
call compute_err1(nO,nV,f_o,f_v,r1,err1)
call compute_err2(nO,nV,f_o,f_v,r2,err2)
tmp_err1 = reshape(err1,(/nO*nV/))
tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
tmp_t1 = reshape(t1 ,(/nO*nV/))
tmp_t2 = reshape(t2 ,(/nO*nO*nV*nV/))
! Update the errors and parameters for the diis
call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t1 ,all_t1 ,nO*nV,cc_diis_depth,nb_iter+1)
call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
call update_all_t (tmp_t2 ,all_t2 ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Gather the different parameters and errors
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,all_err,all_err1,all_err2,cc_diis_depth,&
!$OMP all_t,all_t1,all_t2) &
!$OMP PRIVATE(i,j) &
!$OMP DEFAULT(NONE)
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nV
all_err(i,j) = all_err1(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nO*nV*nV
all_err(i+nO*nV,j) = all_err2(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nV
all_t(i,j) = all_t1(i,j)
enddo
!$OMP END DO NOWAIT
enddo
do j = 1, cc_diis_depth
!$OMP DO
do i = 1, nO*nO*nV*nV
all_t(i+nO*nV,j) = all_t2(i,j)
enddo
!$OMP END DO
enddo
!$OMP END PARALLEL
! Diis
call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t)
! Split the resulting vector
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp_t,tmp_t1,tmp_t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
tmp_t1(i) = tmp_t(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
tmp_t2(i) = tmp_t(i+nO*nV)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Reshape as tensors
t1 = reshape(tmp_t1 ,(/nO,nV/))
t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
! Deallocate
deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2,all_t,all_err)
end
#+end_src
** Update t v3
#+begin_src f90 :comments org :tangle diis.irp.f
subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t)
implicit none
integer, intent(in) :: nO,nV,nb_iter
double precision, intent(in) :: f_o(nO), f_v(nV)
double precision, intent(in) :: r1(nO,nV), r2(nO,nO,nV,nV)
double precision, intent(inout) :: t1(nO*nV), t2(nO*nO*nV*nV)
double precision, intent(inout) :: all_err(nO*nV+nO*nO*nV*nV, cc_diis_depth)
double precision, intent(inout) :: all_t(nO*nV+nO*nO*nV*nV, cc_diis_depth)
double precision, allocatable :: tmp(:)
integer :: i,j
! Allocate
allocate(tmp(nO*nV+nO*nO*nV*nV))
! Compute the errors
call compute_err1(nO,nV,f_o,f_v,r1,tmp(1:nO*nV))
call compute_err2(nO,nV,f_o,f_v,r2,tmp(nO*nV+1:nO*nV+nO*nO*nV*nV))
! Update the errors and parameters for the diis
call update_all_err(tmp,all_err,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp,t1,t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
tmp(i) = t1(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
tmp(i+nO*nV) = t2(i)
enddo
!$OMP END DO
!$OMP END PARALLEL
call update_all_t(tmp,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
! Diis
call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp)
! Split the resulting vector
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tmp,t1,t2) &
!$OMP PRIVATE(i) &
!$OMP DEFAULT(NONE)
!$OMP DO
do i = 1, nO*nV
t1(i) = tmp(i)
enddo
!$OMP END DO NOWAIT
!$OMP DO
do i = 1, nO*nO*nV*nV
t2(i) = tmp(i+nO*nV)
enddo
!$OMP END DO
!$OMP END PARALLEL
! Deallocate
deallocate(tmp)
end
#+end_src

View File

@ -0,0 +1,15 @@
#+begin_src f90 :comments org :tangle energy.irp.f
subroutine det_energy(det,energy)
implicit none
integer(bit_kind), intent(in) :: det
double precision, intent(out) :: energy
call i_H_j(det,det,N_int,energy)
energy = energy + nuclear_repulsion
end
#+end_src

View File

@ -0,0 +1,222 @@
* Guess
** T1
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine guess_t1(nO,nV,f_o,f_v,f_ov,t1)
implicit none
BEGIN_DOC
! Update the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), f_ov(nO,nV)
! inout
double precision, intent(out) :: t1(nO, nV)
! internal
integer :: i,a
if (trim(cc_guess_t1) == 'none') then
t1 = 0d0
else if (trim(cc_guess_t1) == 'MP') then
do a = 1, nV
do i = 1, nO
t1(i,a) = f_ov(i,a) / (f_o(i) - f_v(a) - cc_level_shift_guess)
enddo
enddo
else if (trim(cc_guess_t1) == 'read') then
call read_t1(nO,nV,t1)
else
print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t1)
call abort
endif
end
#+end_src
** T2
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine guess_t2(nO,nV,f_o,f_v,v_oovv,t2)
implicit none
BEGIN_DOC
! Update the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), v_oovv(nO, nO, nV, nV)
! inout
double precision, intent(out) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
if (trim(cc_guess_t2) == 'none') then
t2 = 0d0
else if (trim(cc_guess_t2) == 'MP') then
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
t2(i,j,a,b) = v_oovv(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift_guess)
enddo
enddo
enddo
enddo
else if (trim(cc_guess_t2) == 'read') then
call read_t2(nO,nV,t2)
else
print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t2)
call abort
endif
end
#+end_src
* Write
** T1
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine write_t1(nO,nV,t1)
implicit none
BEGIN_DOC
! Write the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO, nV)
! internal
integer :: i,a
if (cc_write_t1) then
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
do a = 1, nV
do i = 1, nO
write(11,'(F20.12)') t1(i,a)
enddo
enddo
close(11)
endif
end
#+end_src
** T2
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine write_t2(nO,nV,t2)
implicit none
BEGIN_DOC
! Write the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
if (cc_write_t2) then
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
write(11,'(F20.12)') t2(i,j,a,b)
enddo
enddo
enddo
enddo
close(11)
endif
end
#+end_src
* Read
** T1
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine read_t1(nO,nV,t1)
implicit none
BEGIN_DOC
! Read the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(out) :: t1(nO, nV)
! internal
integer :: i,a
logical :: ok
inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
if (.not. ok) then
print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
print*, 'Do a first calculation with cc_write_t1 = True'
print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
call abort
endif
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
do a = 1, nV
do i = 1, nO
read(11,'(F20.12)') t1(i,a)
enddo
enddo
close(11)
end
#+end_src
** T2
#+begin_src f90 :comments org :tangle guess_t.irp.f
subroutine read_t2(nO,nV,t2)
implicit none
BEGIN_DOC
! Read the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(out) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
logical :: ok
inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
if (.not. ok) then
print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
print*, 'Do a first calculation with cc_write_t2 = True'
print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
call abort
endif
open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
read(11,'(F20.12)') t2(i,j,a,b)
enddo
enddo
enddo
enddo
close(11)
end
#+end_src

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,341 @@
* N spin orb
#+begin_src f90 :comments org :tangle occupancy.irp.f
subroutine extract_n_spin(det,n)
implicit none
BEGIN_DOC
! Returns the number of occupied alpha, occupied beta, virtual alpha, virtual beta spin orbitals
! in det without counting the core and deleted orbitals in the format n(nOa,nOb,nVa,nVb)
END_DOC
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: n(4)
integer(bit_kind) :: res(N_int,2)
integer :: i, si
logical :: ok, is_core, is_del
! Init
n = 0
! Loop over the spin
do si = 1, 2
do i = 1, mo_num
call apply_hole(det, si, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
n(si) = n(si) + 1
else
! hole
n(si+2) = n(si+2) + 1
endif
enddo
enddo
!print*,n(1),n(2),n(3),n(4)
end
#+end_src
* List_orb
** Spin
#+begin_src f90 :comments org :tangle occupancy.irp.f
subroutine extract_list_orb_spin(det,nO_m,nV_m,list_occ,list_vir)
implicit none
BEGIN_DOC
! Returns the the list of occupied alpha/beta, virtual alpha/beta spin orbitals
! size(nO_m,1) must be max(nOa,nOb) and size(nV_m,1) must be max(nVa,nVb)
END_DOC
integer, intent(in) :: nO_m, nV_m
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: list_occ(nO_m,2), list_vir(nV_m,2)
integer(bit_kind) :: res(N_int,2)
integer :: i, si, idx_o, idx_v, idx_i, idx_b
logical :: ok, is_core, is_del
list_occ = 0
list_vir = 0
! List of occ/vir alpha/beta
! occ alpha -> list_occ(:,1)
! occ beta -> list_occ(:,2)
! vir alpha -> list_vir(:,1)
! vir beta -> list_vir(:,2)
! Loop over the spin
do si = 1, 2
! tmp idx
idx_o = 1
idx_v = 1
do i = 1, mo_num
call apply_hole(det, si, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
list_occ(idx_o,si) = i
idx_o = idx_o + 1
else
! hole
list_vir(idx_v,si) = i
idx_v = idx_v + 1
endif
enddo
enddo
end
#+end_src
** Space
#+begin_src f90 :comments org :tangle occupancy.irp.f
subroutine extract_list_orb_space(det,nO,nV,list_occ,list_vir)
implicit none
BEGIN_DOC
! Returns the the list of occupied and virtual alpha spin orbitals
END_DOC
integer, intent(in) :: nO, nV
integer(bit_kind), intent(in) :: det(N_int,2)
integer, intent(out) :: list_occ(nO), list_vir(nV)
integer(bit_kind) :: res(N_int,2)
integer :: i, si, idx_o, idx_v, idx_i, idx_b
logical :: ok, is_core, is_del
if (elec_alpha_num /= elec_beta_num) then
print*,'Error elec_alpha_num /= elec_beta_num, impossible to create cc_list_occ and cc_list_vir, abort'
call abort
endif
list_occ = 0
list_vir = 0
! List of occ/vir alpha
! occ alpha -> list_occ(:,1)
! vir alpha -> list_vir(:,1)
! tmp idx
idx_o = 1
idx_v = 1
do i = 1, mo_num
call apply_hole(det, 1, i, res, ok, N_int)
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
if (ok) then
! particle
list_occ(idx_o) = i
idx_o = idx_o + 1
else
! hole
list_vir(idx_v) = i
idx_v = idx_v + 1
endif
enddo
end
#+end_src
** is_core
#+begin_src f90 :comments org :tangle occupancy.irp.f
function is_core(i)
implicit none
BEGIN_DOC
! True if the orbital i is a core orbital
END_DOC
integer, intent(in) :: i
logical :: is_core
integer :: j
! Init
is_core = .False.
! Search
do j = 1, dim_list_core_orb
if (list_core(j) == i) then
is_core = .True.
exit
endif
enddo
end
#+end_src
** is_del
#+begin_src f90 :comments org :tangle occupancy.irp.f
function is_del(i)
implicit none
BEGIN_DOC
! True if the orbital i is a deleted orbital
END_DOC
integer, intent(in) :: i
logical :: is_del
integer :: j
! Init
is_del = .False.
! Search
do j = 1, dim_list_core_orb
if (list_core(j) == i) then
is_del = .True.
exit
endif
enddo
end
#+end_src
* Providers
** N orb
#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
BEGIN_PROVIDER [integer, cc_nO_m]
&BEGIN_PROVIDER [integer, cc_nOa]
&BEGIN_PROVIDER [integer, cc_nOb]
&BEGIN_PROVIDER [integer, cc_nOab]
&BEGIN_PROVIDER [integer, cc_nV_m]
&BEGIN_PROVIDER [integer, cc_nVa]
&BEGIN_PROVIDER [integer, cc_nVb]
&BEGIN_PROVIDER [integer, cc_nVab]
&BEGIN_PROVIDER [integer, cc_n_mo]
&BEGIN_PROVIDER [integer, cc_nO_S, (2)]
&BEGIN_PROVIDER [integer, cc_nV_S, (2)]
implicit none
BEGIN_DOC
! Number of orbitals without core and deleted ones of the cc_ref det in psi_det
! a: alpha, b: beta
! nO_m: max(a,b) occupied
! nOa: nb a occupied
! nOb: nb b occupied
! nOab: nb a+b occupied
! nV_m: max(a,b) virtual
! nVa: nb a virtual
! nVb: nb b virtual
! nVab: nb a+b virtual
END_DOC
integer :: n_spin(4)
! Extract number of occ/vir alpha/beta spin orbitals
call extract_n_spin(psi_det(1,1,cc_ref),n_spin)
cc_nOa = n_spin(1)
cc_nOb = n_spin(2)
cc_nOab = cc_nOa + cc_nOb !n_spin(1) + n_spin(2)
cc_nO_m = max(cc_nOa,cc_nOb) !max(n_spin(1), n_spin(2))
cc_nVa = n_spin(3)
cc_nVb = n_spin(4)
cc_nVab = cc_nVa + cc_nVb !n_spin(3) + n_spin(4)
cc_nV_m = max(cc_nVa,cc_nVb) !max(n_spin(3), n_spin(4))
cc_n_mo = cc_nVa + cc_nVb !n_spin(1) + n_spin(3)
cc_nO_S = (/cc_nOa,cc_nOb/)
cc_nV_S = (/cc_nVa,cc_nVb/)
END_PROVIDER
#+end_src
** List orb
*** General
#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
BEGIN_PROVIDER [integer, cc_list_gen, (cc_n_mo)]
implicit none
BEGIN_DOC
! List of general orbitals without core and deleted ones
END_DOC
integer :: i,j
logical :: is_core, is_del
j = 1
do i = 1, mo_num
! in core ?
if (is_core(i)) cycle
! in del ?
if (is_del(i)) cycle
cc_list_gen(j) = i
j = j+1
enddo
END_PROVIDER
#+end_src
*** Space
#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
BEGIN_PROVIDER [integer, cc_list_occ, (cc_nOa)]
&BEGIN_PROVIDER [integer, cc_list_vir, (cc_nVa)]
implicit none
BEGIN_DOC
! List of occupied and virtual spatial orbitals without core and deleted ones
END_DOC
call extract_list_orb_space(psi_det(1,1,cc_ref),cc_nOa,cc_nVa,cc_list_occ,cc_list_vir)
END_PROVIDER
#+end_src
*** Spin
#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
BEGIN_PROVIDER [integer, cc_list_occ_spin, (cc_nO_m,2)]
&BEGIN_PROVIDER [integer, cc_list_vir_spin, (cc_nV_m,2)]
&BEGIN_PROVIDER [logical, cc_ref_is_open_shell]
implicit none
BEGIN_DOC
! List of occupied and virtual spin orbitals without core and deleted ones
END_DOC
integer :: i
call extract_list_orb_spin(psi_det(1,1,cc_ref),cc_nO_m,cc_nV_m,cc_list_occ_spin,cc_list_vir_spin)
cc_ref_is_open_shell = .False.
do i = 1, cc_nO_m
if (cc_list_occ_spin(i,1) /= cc_list_occ_spin(i,2)) then
cc_ref_is_open_shell = .True.
endif
enddo
END_PROVIDER
#+end_src

178
src/utils_cc/org/phase.org Normal file
View File

@ -0,0 +1,178 @@
#+begin_src f90 :comments org :notangle phase.irp.f
program run
implicit none
integer :: n(2), degree1, degree2, exc(0:2,2,2)
integer, allocatable :: list_anni(:,:), list_crea(:,:)
double precision :: phase1, phase2
integer :: h1,h2,p1,p2,s1,s2,i,j
allocate(list_anni(N_int*bit_kind_size,2))
allocate(list_crea(N_int*bit_kind_size,2))
do i = 1, N_det-1
do j = i+1, N_det
!call print_det(psi_det(1,1,j),N_int)
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree1,phase1,N_int)
call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
!print*,'old',degree1,phase1
!print*,'h1:',h1,'h2:',h2,'s1:',s1,'s2:',s2
!print*,'p1:',p1,'p2:',p2
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree1,N_int)
call get_excitation_general(psi_det(1,1,i),psi_det(1,1,j),degree2,n,list_anni,list_crea,phase2,N_int)
!print*,'new',degree2,phase2
!print*,'ha:',list_anni(1:n(1),1),'hb',list_anni(1:n(2),2)
!print*,'pa:',list_crea(1:n(1),1),'pb',list_crea(1:n(2),2)
!print*,''
if (degree1 /= degree2) then
print*,'Error degree:',degree1,degree2
call abort
endif
if (degree1 <= 2 .and. phase1 /= phase2) then
print*,'Error phase',phase1,phase2
call abort
endif
enddo
enddo
end
#+end_src
** phase
#+begin_src f90 :comments org :tangle phase.irp.f
subroutine get_phase_general(det1,det2,phase,degree,Nint)
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
double precision, intent(out) :: phase
integer, intent(out) :: degree
integer :: n(2)
integer, allocatable :: list_anni(:,:), list_crea(:,:)
allocate(list_anni(N_int*bit_kind_size,2))
allocate(list_crea(N_int*bit_kind_size,2))
call get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
end
#+end_src
** Get excitation general
#+begin_src f90 :comments org :tangle phase.irp.f
subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
double precision, intent(out) :: phase
integer, intent(out) :: list_crea(Nint*bit_kind_size,2)
integer, intent(out) :: list_anni(Nint*bit_kind_size,2)
integer, intent(out) :: degree, n(2)
integer, allocatable :: l1(:,:), l2(:,:)
integer(bit_kind), allocatable :: det_crea(:,:), det_anni(:,:)
integer, allocatable :: pos_anni(:,:), pos_crea(:,:)
integer :: n1(2),n2(2),n_crea(2),n_anni(2),i,j,k,d
allocate(l1(Nint*bit_kind_size,2))
allocate(l2(Nint*bit_kind_size,2))
allocate(det_crea(Nint,2),det_anni(Nint,2))
! 1 111010
! 2 110101
!
!not 1-> 000101
! 2 110101
!and 000101 -> crea
!
! 1 111010
!not 2-> 001010
! 001010 -> anni
do j = 1, 2
do i = 1, Nint
det_crea(i,j) = iand(not(det1(i,j)),det2(i,j))
enddo
enddo
do j = 1, 2
do i = 1, Nint
det_anni(i,j) = iand(det1(i,j),not(det2(i,j)))
enddo
enddo
call bitstring_to_list_ab(det1,l1,n1,Nint)
call bitstring_to_list_ab(det2,l2,n2,Nint)
call bitstring_to_list_ab(det_crea,list_crea,n_crea,Nint)
call bitstring_to_list_ab(det_anni,list_anni,n_anni,Nint)
do i = 1, 2
if (n_crea(i) /= n_anni(i)) then
print*,'Well, it seems we have a problem here...'
call abort
endif
enddo
!1 11110011001 1 2 3 4 7 8 11
!pos 1 2 3 4 5 6 7
!2 11100101011 1 2 3 6 8 10 11
!anni 00010010000 4 7
!pos 4 5
!crea 00000100010 6 10
!pos 4 6
!4 -> 6 pos(4 -> 4)
!7 -> 10 pos(5 -> 6)
n = n_anni
degree = n_anni(1) + n_anni(2)
allocate(pos_anni(max(n(1),n(2)),2))
allocate(pos_crea(max(n(1),n(2)),2))
! Search pos anni
do j = 1, 2
k = 1
do i = 1, n1(j)
if (l1(i,j) /= list_anni(k,j)) cycle
pos_anni(k,j) = i
k = k + 1
enddo
enddo
! Search pos crea
do j = 1, 2
k = 1
do i = 1, n2(j)
if (l2(i,j) /= list_crea(k,j)) cycle
pos_crea(k,j) = i
k = k + 1
enddo
enddo
! Distance between the ith anni and the ith crea op
! By doing so there is no crossing between the different pairs of anni/crea
! and the phase is determined by the sum of the distances
! -> (-1)^{sum of the distances}
d = 0
do j = 1, 2
do i = 1, n(j)
d = d + abs(pos_anni(i,j) - pos_crea(i,j))
enddo
enddo
phase = dble((-1)**d)
! Debug
!print*,l2(1:n2(1),1)
!print*,l2(1:n2(2),2)
!!call print_det(det1,Nint)
!!call print_det(det2,Nint)
!print*,phase
!print*,''
end
#+end_src

View File

@ -0,0 +1,33 @@
#+begin_src f90 :comments org :tangle print_wf_qp_edit.irp.f
program run
implicit none
read_wf = .true.
touch read_wf
call print_wf_qp_edit()
end
#+end_src
#+begin_src f90 :comments org :tangle print_wf_qp_edit.irp.f
subroutine print_wf_qp_edit()
implicit none
BEGIN_DOC
! Print the psi_det wave function up to n_det_qp_edit
END_DOC
integer :: i
do i = 1, n_det_qp_edit
print*,i
write(*,'(100(1pE12.4))') psi_coef(i,:)
call print_det(psi_det(1,1,i),N_int)
print*,''
enddo
end
#+end_src

View File

@ -0,0 +1,76 @@
* T1
#+begin_src f90 :comments org :tangle update_t.irp.f
subroutine update_t1(nO,nV,f_o,f_v,r1,t1)
implicit none
BEGIN_DOC
! Update the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r1(nO, nV)
! inout
double precision, intent(inout) :: t1(nO, nV)
! internal
integer :: i,a
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,t1,r1,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,a) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(1)
do a = 1, nV
do i = 1, nO
t1(i,a) = t1(i,a) - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
#+end_src
* T2
#+begin_src f90 :comments org :tangle update_t.irp.f
subroutine update_t2(nO,nV,f_o,f_v,r2,t2)
implicit none
BEGIN_DOC
! Update the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r2(nO, nO, nV, nV)
! inout
double precision, intent(inout) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,t2,r2,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
t2(i,j,a,b) = t2(i,j,a,b) - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
#+end_src

135
src/utils_cc/phase.irp.f Normal file
View File

@ -0,0 +1,135 @@
! phase
subroutine get_phase_general(det1,det2,phase,degree,Nint)
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
double precision, intent(out) :: phase
integer, intent(out) :: degree
integer :: n(2)
integer, allocatable :: list_anni(:,:), list_crea(:,:)
allocate(list_anni(N_int*bit_kind_size,2))
allocate(list_crea(N_int*bit_kind_size,2))
call get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
end
! Get excitation general
subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
double precision, intent(out) :: phase
integer, intent(out) :: list_crea(Nint*bit_kind_size,2)
integer, intent(out) :: list_anni(Nint*bit_kind_size,2)
integer, intent(out) :: degree, n(2)
integer, allocatable :: l1(:,:), l2(:,:)
integer(bit_kind), allocatable :: det_crea(:,:), det_anni(:,:)
integer, allocatable :: pos_anni(:,:), pos_crea(:,:)
integer :: n1(2),n2(2),n_crea(2),n_anni(2),i,j,k,d
allocate(l1(Nint*bit_kind_size,2))
allocate(l2(Nint*bit_kind_size,2))
allocate(det_crea(Nint,2),det_anni(Nint,2))
! 1 111010
! 2 110101
!
!not 1-> 000101
! 2 110101
!and 000101 -> crea
!
! 1 111010
!not 2-> 001010
! 001010 -> anni
do j = 1, 2
do i = 1, Nint
det_crea(i,j) = iand(not(det1(i,j)),det2(i,j))
enddo
enddo
do j = 1, 2
do i = 1, Nint
det_anni(i,j) = iand(det1(i,j),not(det2(i,j)))
enddo
enddo
call bitstring_to_list_ab(det1,l1,n1,Nint)
call bitstring_to_list_ab(det2,l2,n2,Nint)
call bitstring_to_list_ab(det_crea,list_crea,n_crea,Nint)
call bitstring_to_list_ab(det_anni,list_anni,n_anni,Nint)
do i = 1, 2
if (n_crea(i) /= n_anni(i)) then
print*,'Well, it seems we have a problem here...'
call abort
endif
enddo
!1 11110011001 1 2 3 4 7 8 11
!pos 1 2 3 4 5 6 7
!2 11100101011 1 2 3 6 8 10 11
!anni 00010010000 4 7
!pos 4 5
!crea 00000100010 6 10
!pos 4 6
!4 -> 6 pos(4 -> 4)
!7 -> 10 pos(5 -> 6)
n = n_anni
degree = n_anni(1) + n_anni(2)
allocate(pos_anni(max(n(1),n(2)),2))
allocate(pos_crea(max(n(1),n(2)),2))
! Search pos anni
do j = 1, 2
k = 1
do i = 1, n1(j)
if (l1(i,j) /= list_anni(k,j)) cycle
pos_anni(k,j) = i
k = k + 1
enddo
enddo
! Search pos crea
do j = 1, 2
k = 1
do i = 1, n2(j)
if (l2(i,j) /= list_crea(k,j)) cycle
pos_crea(k,j) = i
k = k + 1
enddo
enddo
! Distance between the ith anni and the ith crea op
! By doing so there is no crossing between the different pairs of anni/crea
! and the phase is determined by the sum of the distances
! -> (-1)^{sum of the distances}
d = 0
do j = 1, 2
do i = 1, n(j)
d = d + abs(pos_anni(i,j) - pos_crea(i,j))
enddo
enddo
phase = dble((-1)**d)
! Debug
!print*,l2(1:n2(1),1)
!print*,l2(1:n2(2),2)
!!call print_det(det1,Nint)
!!call print_det(det2,Nint)
!print*,phase
!print*,''
end

View File

@ -0,0 +1,29 @@
program run
implicit none
read_wf = .true.
touch read_wf
call print_wf_qp_edit()
end
subroutine print_wf_qp_edit()
implicit none
BEGIN_DOC
! Print the psi_det wave function up to n_det_qp_edit
END_DOC
integer :: i
do i = 1, n_det_qp_edit
print*,i
write(*,'(100(1pE12.4))') psi_coef(i,:)
call print_det(psi_det(1,1,i),N_int)
print*,''
enddo
end

View File

@ -0,0 +1,73 @@
! T1
subroutine update_t1(nO,nV,f_o,f_v,r1,t1)
implicit none
BEGIN_DOC
! Update the T1 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r1(nO, nV)
! inout
double precision, intent(inout) :: t1(nO, nV)
! internal
integer :: i,a
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,t1,r1,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,a) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(1)
do a = 1, nV
do i = 1, nO
t1(i,a) = t1(i,a) - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
! T2
subroutine update_t2(nO,nV,f_o,f_v,r2,t2)
implicit none
BEGIN_DOC
! Update the T2 amplitudes for CC
END_DOC
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: f_o(nO), f_v(nV), r2(nO, nO, nV, nV)
! inout
double precision, intent(inout) :: t2(nO, nO, nV, nV)
! internal
integer :: i,j,a,b
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,t2,r2,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
t2(i,j,a,b) = t2(i,j,a,b) - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end