4
1
mirror of https://github.com/pfloos/quack synced 2024-11-06 22:24:03 +01:00

spin flip for CIS and TDHF almost done

This commit is contained in:
Pierre-Francois Loos 2020-10-05 16:58:19 +02:00
parent 75ece830a6
commit 3d2ca7529f
46 changed files with 965 additions and 192 deletions

30
examples/basis.B.cc-pvdz Normal file
View File

@ -0,0 +1,30 @@
1 6
S 8
1 4570.0000000 0.0006960
2 685.9000000 0.0053530
3 156.5000000 0.0271340
4 44.4700000 0.1013800
5 14.4800000 0.2720550
6 5.1310000 0.4484030
7 1.8980000 0.2901230
8 0.3329000 0.0143220
S 8
1 4570.0000000 -0.0001390
2 685.9000000 -0.0010970
3 156.5000000 -0.0054440
4 44.4700000 -0.0219160
5 14.4800000 -0.0597510
6 5.1310000 -0.1387320
7 1.8980000 -0.1314820
8 0.3329000 0.5395260
S 1
1 0.1043000 1.0000000
P 3
1 6.0010000 0.0354810
2 1.2410000 0.1980720
3 0.3364000 0.5052300
P 1
1 0.0953800 1.0000000
D 1
1 0.3430000 1.0000000

22
examples/basis.Be.6-31g Normal file
View File

@ -0,0 +1,22 @@
1 5
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3
1 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000

View File

@ -0,0 +1,23 @@
1 6
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3
1 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
D 1
1 0.4000000 1.0000000

32
examples/basis.C.cc-pvdz Normal file
View File

@ -0,0 +1,32 @@
1 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
S 1
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
P 1
1 1.517000E-01 1.000000E+00
D 1
1 5.500000E-01 1.0000000

64
examples/basis.C2.cc-pvdz Normal file
View File

@ -0,0 +1,64 @@
1 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
S 1
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
P 1
1 1.517000E-01 1.000000E+00
D 1
1 5.500000E-01 1.0000000
2 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
S 1
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
P 1
1 1.517000E-01 1.000000E+00
D 1
1 5.500000E-01 1.0000000

View File

@ -0,0 +1,49 @@
1 9
S 8
1 6665.0000000 0.0006920
2 1000.0000000 0.0053290
3 228.0000000 0.0270770
4 64.7100000 0.1017180
5 21.0600000 0.2747400
6 7.4950000 0.4485640
7 2.7970000 0.2850740
8 0.5215000 0.0152040
S 8
1 6665.0000000 -0.0001460
2 1000.0000000 -0.0011540
3 228.0000000 -0.0057250
4 64.7100000 -0.0233120
5 21.0600000 -0.0639550
6 7.4950000 -0.1499810
7 2.7970000 -0.1272620
8 0.5215000 0.5445290
S 1
1 0.1596000 1.0000000
S 1
1 0.0469000 1.0000000
P 3
1 9.4390000 0.0381090
2 2.0020000 0.2094800
3 0.5456000 0.5085570
P 1
1 0.1517000 1.0000000
P 1
1 0.0404100 1.0000000
D 1
1 0.5500000 1.0000000
D 1
1 0.1510000 1.0000000
2 5
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
S 1
1 0.0297400 1.0000000
P 1
1 0.7270000 1.0000000
P 1
1 0.1410000 1.0000000

38
examples/basis.CH.cc-pvdz Normal file
View File

@ -0,0 +1,38 @@
1 6
S 8
1 6665.0000000 0.0006920
2 1000.0000000 0.0053290
3 228.0000000 0.0270770
4 64.7100000 0.1017180
5 21.0600000 0.2747400
6 7.4950000 0.4485640
7 2.7970000 0.2850740
8 0.5215000 0.0152040
S 8
1 6665.0000000 -0.0001460
2 1000.0000000 -0.0011540
3 228.0000000 -0.0057250
4 64.7100000 -0.0233120
5 21.0600000 -0.0639550
6 7.4950000 -0.1499810
7 2.7970000 -0.1272620
8 0.5215000 0.5445290
S 1
1 0.1596000 1.0000000
P 3
1 9.4390000 0.0381090
2 2.0020000 0.2094800
3 0.5456000 0.5085570
P 1
1 0.1517000 1.0000000
D 1
1 0.5500000 1.0000000
2 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000

28
examples/basis.CO.sto-3g Normal file
View File

@ -0,0 +1,28 @@
1 3
S 3
1 71.6168370 0.15432897
2 13.0450960 0.53532814
3 3.5305122 0.44463454
S 3
1 2.9412494 -0.0999672
2 0.6834831 0.3995128
3 0.2222899 0.7001155
P 3
1 2.9412494 0.1559163
2 0.6834831 0.6076837
3 0.2222899 0.3919574
2 3
S 3
1 130.7093200 0.15432897
2 23.8088610 0.53532814
3 6.4436083 0.44463454
S 3
1 5.0331513 -0.0999672
2 1.1695961 0.3995128
3 0.3803890 0.7001155
P 3
1 5.0331513 0.1559163
2 1.1695961 0.6076837
3 0.3803890 0.3919574

29
examples/basis.F.cc-pvdz Normal file
View File

@ -0,0 +1,29 @@
1 6
S 8
1 14710.0000000 0.0007210
2 2207.0000000 0.0055530
3 502.8000000 0.0282670
4 142.6000000 0.1064440
5 46.4700000 0.2868140
6 16.7000000 0.4486410
7 6.3560000 0.2647610
8 1.3160000 0.0153330
S 8
1 14710.0000000 -0.0001650
2 2207.0000000 -0.0013080
3 502.8000000 -0.0064950
4 142.6000000 -0.0266910
5 46.4700000 -0.0736900
6 16.7000000 -0.1707760
7 6.3560000 -0.1123270
8 1.3160000 0.5628140
S 1
1 0.3897000 1.0000000
P 3
1 22.6700000 0.0448780
2 4.9770000 0.2357180
3 1.3470000 0.5085210
P 1
1 0.3471000 1.0000000
D 1
1 1.6400000 1.0000000

9
examples/basis.H.cc-pvdz Normal file
View File

@ -0,0 +1,9 @@
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000

5
examples/basis.H.sto-3g Normal file
View File

@ -0,0 +1,5 @@
1 1
S 3
1 3.42525091 0.15432897
2 0.62391373 0.53532814
3 0.16885540 0.44463454

18
examples/basis.H2.6-311g Normal file
View File

@ -0,0 +1,18 @@
1 3
S 3
1 33.8650000 0.0254938
2 5.0947900 0.1903730
3 1.1587900 0.8521610
S 1
1 0.3258400 1.0000000
S 1
1 0.1027410 1.0000000
2 3
S 3
1 33.8650000 0.0254938
2 5.0947900 0.1903730
3 1.1587900 0.8521610
S 1
1 0.3258400 1.0000000
S 1
1 0.1027410 1.0000000

6
examples/basis.H2.pipis Normal file
View File

@ -0,0 +1,6 @@
1 1
P 1
1 0.7270000 1.0000000
2 1
P 1
1 0.7270000 1.0000000

View File

@ -0,0 +1,11 @@
1 1
S 3
1 6.36242139 0.15432897
2 1.15892300 0.53532814
3 0.31364979 0.44463454
2 1
S 3
1 3.42525091 0.15432897
2 0.62391373 0.53532814
3 0.16885540 0.44463454

View File

@ -0,0 +1,30 @@
1 6
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
S 1
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
P 1
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000

30
examples/basis.Li.cc-pvdz Normal file
View File

@ -0,0 +1,30 @@
1 6
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
S 1
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
P 1
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000

14
examples/basis.Li.sto-3g Normal file
View File

@ -0,0 +1,14 @@
1 3
S 3
1 16.1195750 0.15432897
2 2.9362007 0.53532814
3 0.7946505 0.44463454
S 3
1 0.6362897 -0.0999672
2 0.1478601 0.3995128
3 0.0480887 0.7001155
P 3
1 0.6362897 0.1559163
2 0.1478601 0.6076837
3 0.0480887 0.3919574

21
examples/basis.Ne.6-31g Normal file
View File

@ -0,0 +1,21 @@
1 5
S 6
1 8425.8515300 0.0018843481
2 1268.5194000 0.0143368994
3 289.6214140 0.0701096233
4 81.8590040 0.2373732660
5 26.2515079 0.4730071260
6 9.09472051 0.3484012410
S 3
1 26.5321310 -0.1071183
2 6.1017550 -0.1461638
3 1.6962715 1.1277735
P 3
1 26.5321310 0.0719096
2 6.1017550 0.3495134
3 1.6962715 0.7199405
S 1
1 0.4458187 1.0000000
P 1
1 0.4458187 1.0000000

View File

@ -0,0 +1,23 @@
1 6
S 6
1 8425.8515300 0.0018843481
2 1268.5194000 0.0143368994
3 289.6214140 0.0701096233
4 81.8590040 0.2373732660
5 26.2515079 0.4730071260
6 9.09472051 0.3484012410
S 3
1 26.5321310 -0.1071183
2 6.1017550 -0.1461638
3 1.6962715 1.1277735
P 3
1 26.5321310 0.0719096
2 6.1017550 0.3495134
3 1.6962715 0.7199405
S 1
1 0.4458187 1.0000000
P 1
1 0.4458187 1.0000000
D 1
1 0.8000000 1.0000000

4
examples/molecule.B Normal file
View File

@ -0,0 +1,4 @@
# nAt nEla nElb nCore nRyd
1 3 2 0 0
# Znuc x y z
B 0. 0. 0.

4
examples/molecule.C Normal file
View File

@ -0,0 +1,4 @@
# nAt nEla nElb nCore nRyd
1 5 1 0 0
# Znuc x y z
C 0. 0. 0.

5
examples/molecule.CH Normal file
View File

@ -0,0 +1,5 @@
# nAt nEla nElb nCore nRyd
2 4 3 0 0
# Znuc x y z
C 0. 0. -0.16245872
H 0. 0. 1.93436816

4
examples/molecule.F Normal file
View File

@ -0,0 +1,4 @@
# nAt nEla nElb nCore nRyd
1 5 4 0 0
# Znuc x y z
F 0. 0. 0.

4
examples/molecule.H Normal file
View File

@ -0,0 +1,4 @@
# nAt nEla nElb nCore nRyd
1 1 0 0 0
# Znuc x y z
H 0. 0. 0.

4
examples/molecule.Li+ Normal file
View File

@ -0,0 +1,4 @@
# nAt nEla nElb nCore nRyd
1 1 1 0 0
# Znuc x y z
Li 0.0 0.0 0.0

View File

@ -1,5 +1,5 @@
# RHF UHF MOM
T F F
F T F
# MP2* MP3 MP2-F12
F F F
# CCD CCSD CCSD(T)
@ -7,9 +7,9 @@
# drCCD rCCD lCCD pCCD
F F F F
# CIS* CIS(D) CID CISD
F F F F
T F F F
# RPA* RPAx* ppRPA
F T F
F F F
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0* evGW* qsGW

View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: singlet triplet spin_conserved spin_flip TDA
T T T F T
T T T T F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0

View File

@ -1,5 +1,5 @@
subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
dipole_int_aa,dipole_int_bb,eHF)
dipole_int_aa,dipole_int_bb,eHF,cHF,S)
! Perform configuration interaction single calculation`
@ -17,6 +17,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -80,8 +82,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
call diagonalize_matrix(nS_sc,A_sc,Omega_sc)
call print_excitation('UCIS ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,transpose(A_sc),transpose(A_sc))
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
cHF,S,Omega_sc,transpose(A_sc),transpose(A_sc))
if(dump_trans) then
print*,'Spin-conserved CIS transition vectors'
@ -118,14 +120,14 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
write(*,*)
endif
! call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
allocate(order(nS_sf))
call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf)
call quick_sort(Omega_sf,order(:),nS_sf)
call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
! allocate(order(nS_sf))
! call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf)
! call quick_sort(Omega_sf,order(:),nS_sf)
call print_excitation('UCIS ',6,nS_sf,Omega_sf)
call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
Omega_sf,transpose(A_sf),transpose(A_sf))
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
cHF,S,Omega_sf,transpose(A_sf),transpose(A_sf))
if(dump_trans) then
print*,'Spin-flip CIS transition vectors'

View File

@ -1,4 +1,4 @@
subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
! Print one- and two-electron energies and other stuff for UHF calculation
@ -7,7 +7,7 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: ENuc
@ -19,14 +19,12 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
double precision,intent(in) :: dipole(ncart)
integer :: ixyz
integer :: i,j
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
double precision :: S2_exact
double precision :: S2
integer :: spin_state
double precision :: S_exact,S2_exact
double precision :: S,S2
! HOMO and LUMO
@ -47,9 +45,10 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
end do
S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2)
S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
spin_state = nO(1) - nO(2) + 1
S_exact = 0.5d0*dble(nO(1) - nO(2))
S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
! Dump results
@ -91,7 +90,8 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,I6)') ' 2S+1 :',spin_state
write(*,'(A40,F13.6)') ' S (exact) :',2d0*S_exact + 1d0
write(*,'(A40,F13.6)') ' S :',2d0*S + 1d0
write(*,'(A40,F13.6)') ' <S**2> (exact) :',S2_exact
write(*,'(A40,F13.6)') ' <S**2> :',S2
write(*,'(A60)') '-------------------------------------------------'

View File

@ -69,9 +69,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
do ia=1,nS
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq)
call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv)
@ -83,9 +83,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: negative excitations!!')
do ia=1,nS
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
Omega = sqrt(Omega)

View File

@ -1,4 +1,4 @@
subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os)
! Compute linear response
@ -13,6 +13,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: maxS
double precision :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: XpY(nS,nS)
@ -20,7 +21,6 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
@ -32,7 +32,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
! Memory allocation
allocate(f(nS,ncart))
allocate(f(maxS,ncart))
! Initialization
@ -40,7 +40,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
! Compute dipole moments and oscillator strengths
do ia=1,nS
do ia=1,maxS
do ixyz=1,ncart
jb = 0
do j=nC+1,nO
@ -53,24 +53,19 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
end do
f(:,:) = sqrt(2d0)*f(:,:)
do ia=1,nS
do ia=1,maxS
os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
end do
if(debug) then
write(*,*) '------------------------'
write(*,*) ' Dipole moments (X Y Z) '
write(*,*) '------------------------'
call matout(nS,ncart,f)
write(*,*) '---------------------------------------------------------------'
write(*,*) ' Transition dipole moment (au) '
write(*,*) '---------------------------------------------------------------'
write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.'
write(*,*) '---------------------------------------------------------------'
do ia=1,maxS
write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia)
end do
write(*,*) '---------------------------------------------------------------'
write(*,*)
write(*,*) '----------------------'
write(*,*) ' Oscillator strengths '
write(*,*) '----------------------'
call matout(nS,1,os)
write(*,*)
end if
end subroutine oscillator_strength

View File

@ -21,10 +21,8 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
integer,parameter :: maxS = 10
integer :: ia,jb,j,b
integer :: maxS = 10
double precision :: S2
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:)
@ -33,20 +31,29 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,
! Memory allocation
allocate(X(nS),Y(nS),os(nS))
maxS = min(nS,maxS)
allocate(X(nS),Y(nS),os(maxS))
! Compute oscillator strengths
os(:) = 0d0
if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os)
if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os)
! Print details about excitations
do ia=1,min(nS,maxS)
do ia=1,maxS
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
! <S**2> values
if(spin_allowed) then
S2 = 0d0
else
S2 = 2d0
end if
print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2

View File

@ -1,5 +1,5 @@
subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, &
Omega,XpY,XmY)
subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, &
c,S,Omega,XpY,XmY)
! Print transition vectors for linear response calculation
@ -8,7 +8,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
! Input variables
logical,intent(in) :: spin_allowed
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
@ -20,54 +20,50 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
integer,intent(in) :: nSt
double precision :: dipole_int_aa(nBas,nBas,ncart)
double precision :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt)
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
integer :: ispin
integer,parameter :: maxS = 10
double precision :: S2
integer :: ia,jb,j,b
integer :: maxS = 10
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:)
double precision,allocatable :: Y(:)
double precision,allocatable :: f(:,:)
double precision,allocatable :: os(:)
double precision,allocatable :: S2(:)
! Memory allocation
allocate(X(nSt),Y(nSt),os(nSt))
maxS = min(nSt,maxS)
allocate(X(nSt),Y(nSt),os(maxS),S2(maxS))
! Compute oscillator strengths
os(:) = 0d0
if(spin_allowed) call unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, &
if(ispin == 1) call unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, &
dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os)
! Print details about excitations
! Compute <S**2>
do ia=1,min(nSt,maxS)
call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
! Print details about spin-conserved excitations
if(ispin == 1) then
do ia=1,maxS
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
S2 = (nO(1) - nO(2))/2d0
S2 = 2d0*S2+1d0
S2 = 0.0d0
do jb=1,nSa
S2 = S2 + 4d0*(X(jb)**2 + Y(jb)**2)
end do
do jb=1,nSb
S2 = S2 - 4d0*(X(nSa+jb)**2 + Y(nSa+jb)**2)
end do
print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia)
print*,'-------------------------------------------------------------'
! Spin-up transitions
@ -109,6 +105,64 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
end do
end if
! Print details about spin-flip excitations
if(ispin == 2) then
do ia=1,maxS
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia)
print*,'-------------------------------------------------------------'
! Spin-up transitions
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(jb)
end do
end do
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(jb)
end do
end do
! Spin-down transitions
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(nSa+jb)
end do
end do
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(nSa+jb)
end do
end do
write(*,*)
end do
end if
! Thomas-Reiche-Kuhn sum rule
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))

View File

@ -0,0 +1,177 @@
subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
! Compute <S**2> for linear response excited states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: maxS
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt)
! Local variables
integer :: m
integer :: ia,i,a
double precision :: S2_exact
double precision :: S2_gs
double precision,allocatable :: Xa(:,:), Xb(:,:), Ya(:,:), Yb(:,:)
double precision,allocatable :: Xat(:,:),Xbt(:,:),Yat(:,:),Ybt(:,:)
double precision,allocatable :: OO(:,:), OV(:,:), VO(:,:), VV(:,:)
double precision,allocatable :: OOt(:,:),OVt(:,:),VOt(:,:),VVt(:,:)
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: S2(maxS)
! Memory allocation
allocate(OO(nO(1)-nC(1),nO(2)-nC(2)), OV(nO(1)-nC(1),nV(2)-nR(2)), VO(nV(1)-nR(1),nO(2)-nC(2)), VV(nV(1)-nR(1),nV(2)-nR(2)), &
OOt(nO(2)-nC(2),nO(1)-nC(1)),OVt(nV(2)-nR(2),nO(1)-nC(1)),VOt(nO(2)-nC(2),nV(1)-nR(1)),VVt(nV(2)-nR(2),nV(1)-nR(1)))
! Overlap matrix between spin-up and spin-down orbitals
OO(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nC(2)+1:nO(2) ,2)))
OV(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2)))
VO(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nC(2)+1:nO(2) ,2)))
VV(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2)))
OOt(:,:) = transpose(OO(:,:))
OVt(:,:) = transpose(OV(:,:))
VOt(:,:) = transpose(VO(:,:))
VVt(:,:) = transpose(VV(:,:))
!-------------------------!
! <S**2> for ground state !
!-------------------------!
S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
S2_gs = S2_exact + nO(2) - sum(OO(:,:)**2)
!------------------------------------------!
! <S**2> for spin-conserved-excited states !
!------------------------------------------!
if(ispin == 1) then
allocate(Xa(nO(1)-nC(1),nV(1)-nR(1)), Ya(nO(1)-nC(1),nV(1)-nR(1)), Xb(nO(2)-nC(2),nV(2)-nR(2)), Yb(nO(2)-nC(2),nV(2)-nR(2)), &
Xat(nV(1)-nR(1),nO(1)-nC(1)),Yat(nV(1)-nR(1),nO(1)-nC(1)),Xbt(nV(2)-nR(2),nO(2)-nC(2)),Ybt(nV(2)-nR(2),nO(2)-nC(2)))
do m=1,maxS
ia = 0
do i=nC(1)+1,nO(1)
do a=1,nV(1)-nR(1)
ia = ia + 1
Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia))
Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia))
end do
end do
ia = 0
do i=nC(2)+1,nO(2)
do a=1,nV(2)-nR(2)
ia = ia + 1
Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia))
Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia))
end do
end do
Xat(:,:) = transpose(Xa(:,:))
Xbt(:,:) = transpose(Xb(:,:))
Yat(:,:) = transpose(Ya(:,:))
Ybt(:,:) = transpose(Yb(:,:))
S2(m) = S2_gs &
+ trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) &
+ trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
- trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) &
- trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) &
- 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) &
- 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) &
- 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) &
- trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) &
- trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) &
+ trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) &
+ trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) &
+ 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt))))
end do
end if
!------------------------------------------!
! <S**2> for spin-conserved-excited states !
!------------------------------------------!
if(ispin == 2) then
allocate(Xa(nO(1)-nC(1),nV(2)-nR(2)), Ya(nO(1)-nC(1),nV(2)-nR(2)), Xb(nO(2)-nC(2),nV(1)-nR(1)), Yb(nO(2)-nC(2),nV(1)-nR(1)), &
Xat(nV(2)-nR(2),nO(1)-nC(1)),Yat(nV(2)-nR(2),nO(1)-nC(1)),Xbt(nV(1)-nR(1),nO(2)-nC(2)),Ybt(nV(1)-nR(1),nO(2)-nC(2)))
do m=1,maxS
ia = 0
do i=nC(1)+1,nO(1)
do a=1,nV(2)-nR(2)
ia = ia + 1
Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia))
Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia))
end do
end do
ia = 0
do i=nC(2)+1,nO(2)
do a=1,nV(1)-nR(1)
ia = ia + 1
Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia))
Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia))
end do
end do
Xat(:,:) = transpose(Xa(:,:))
Xbt(:,:) = transpose(Xb(:,:))
Yat(:,:) = transpose(Ya(:,:))
Ybt(:,:) = transpose(Yb(:,:))
S2(m) = 0d0
! S2(m) = S2_gs &
! + trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) &
! + trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
! - trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) &
! - trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) &
! - 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) &
!
! - 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) &
! - 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) &
!
! - trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) &
! - trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) &
! + trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) &
! + trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) &
! + 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt))))
end do
end if
end subroutine unrestricted_S2_expval

View File

@ -35,7 +35,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
! Local variables
integer :: ia
double precision :: trace_matrix
double precision,external :: trace_matrix
double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:)
double precision,allocatable :: ApB(:,:)
@ -90,9 +90,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
do ia=1,nSt
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
! do ia=1,nSt
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq)
call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv)
@ -104,9 +104,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: negative excitations!!')
do ia=1,nSt
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
! do ia=1,nSt
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
Omega = sqrt(Omega)

View File

@ -121,9 +121,6 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
end do
end do
print*,'nSa,nSb,nSt',nSa,nSb,nSt
call matout(nSt,nSt,A_lr)
end if
!-----------------------------------------------
@ -144,10 +141,8 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)'
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i)
- (1d0 - delta_dRPA)*lambda*ERI_aabb(i,b,j,a)
end do
end do
end do
@ -165,17 +160,13 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb)
- (1d0 - delta_dRPA)*lambda*ERI_aabb(b,i,a,j)
end do
end do
end do
end do
print*,'nSa,nSb,nSt',nSa,nSb,nSt
call matout(nSt,nSt,A_lr)
end if

View File

@ -139,7 +139,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j)
B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,j,b,a)
end do
end do
@ -157,7 +157,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a)
B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(j,i,a,b)
end do
end do

View File

@ -1,4 +1,4 @@
subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os)
subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os)
! Compute linear response
@ -17,6 +17,7 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: maxS
double precision :: dipole_int_aa(nBas,nBas,ncart)
double precision :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: Omega(nSt)
@ -25,7 +26,6 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
@ -33,11 +33,11 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo
! Output variables
double precision :: os(nSt)
double precision :: os(maxS)
! Memory allocation
allocate(f(nSt,ncart))
allocate(f(maxS,ncart))
! Initialization
@ -45,7 +45,7 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo
! Compute dipole moments and oscillator strengths
do ia=1,nSt
do ia=1,maxS
do ixyz=1,ncart
jb = 0
@ -67,24 +67,19 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo
end do
end do
do ia=1,nSt
do ia=1,maxS
os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
end do
if(debug) then
write(*,*) '------------------------'
write(*,*) ' Dipole moments (X Y Z) '
write(*,*) '------------------------'
call matout(nS,ncart,f)
write(*,*) '---------------------------------------------------------------'
write(*,*) ' Transition dipole moment (au) '
write(*,*) '---------------------------------------------------------------'
write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.'
write(*,*) '---------------------------------------------------------------'
do ia=1,maxS
write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia)
end do
write(*,*) '---------------------------------------------------------------'
write(*,*)
write(*,*) '----------------------'
write(*,*) ' Oscillator strengths '
write(*,*) '----------------------'
call matout(nS,1,os)
write(*,*)
end if
end subroutine unrestricted_oscillator_strength

View File

@ -98,8 +98,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
! call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
! OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
@ -131,6 +131,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
! Compute spin-flip BSE excitation energies
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)

View File

@ -117,7 +117,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi
end do
end do
@ -141,7 +141,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi
end do
end do

View File

@ -118,7 +118,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi
end do
end do
@ -142,7 +142,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi
end do
end do

View File

@ -52,7 +52,7 @@ debug:
.DEFAULT_GOAL := default
.PHONY: default debug
.PRECIOUS: $(LDIR)/%.a %/Makefile
.PRECIOUS: $(wildcard $(LDIR)/*.a) $(wildcard $(SDIR)/*/Makefile)
clean:
rm -f -- $(patsubst %/,$(LDIR)/%.a,$(wildcard */)) ; \

View File

@ -356,7 +356,8 @@ program QuAcK
! Memory allocation
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas), &
ERI_MO_bbbb(nBas,nBas,nBas,nBas),ERI_MO_abab(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
@ -382,10 +383,6 @@ program QuAcK
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
if(spin_flip) then
allocate(ERI_MO_abab(nBas,nBas,nBas,nBas))
! 4-index transform for (ab|ab) block
bra1 = 1
@ -394,9 +391,6 @@ program QuAcK
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab)
end if
else
! Memory allocation
@ -618,7 +612,7 @@ program QuAcK
if(unrestricted) then
call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb, &
ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S)
else
@ -675,7 +669,7 @@ program QuAcK
if(unrestricted) then
call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S)
else
@ -700,7 +694,7 @@ program QuAcK
if(unrestricted) then
call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S)
else

View File

@ -1,5 +1,5 @@
subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e,c,S)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@ -24,6 +24,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -87,8 +89,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPAx ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,XpY_sc,XmY_sc)
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sc,XpY_sc,XmY_sc)
deallocate(Omega_sc,XpY_sc,XmY_sc)
@ -111,8 +113,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPAx ',6,nS_sf,Omega_sf)
call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
Omega_sf,XpY_sf,XmY_sf)
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sf,XpY_sf,XmY_sf)
deallocate(Omega_sf,XpY_sf,XmY_sf)

View File

@ -1,5 +1,5 @@
subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e,c,S)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@ -24,6 +24,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -86,8 +88,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPA ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,XpY_sc,XmY_sc)
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sc,XpY_sc,XmY_sc)
endif
@ -109,8 +111,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPA ',6,nS_sf,Omega_sf)
call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
Omega_sf,XpY_sf,XmY_sf)
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sf,XpY_sf,XmY_sf)
endif

View File

@ -0,0 +1,45 @@
subroutine print_excitation(method,ispin,nS,Omega)
! Print excitation energies for a given spin manifold
implicit none
include 'parameters.h'
! Input variables
character*12,intent(in) :: method
integer,intent(in) :: ispin,nS
double precision,intent(in) :: Omega(nS)
! Local variables
character*14 :: spin_manifold
integer,parameter :: maxS = 32
integer :: ia
if(ispin == 1) spin_manifold = 'singlet'
if(ispin == 2) spin_manifold = 'triplet'
if(ispin == 3) spin_manifold = 'alpha-beta'
if(ispin == 4) spin_manifold = 'alpha-alpha'
if(ispin == 5) spin_manifold = 'spin-conserved'
if(ispin == 6) spin_manifold = 'spin-flip'
write(*,*)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
write(*,*)'-------------------------------------------------------------'
do ia=1,min(nS,maxS)
write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
'|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_excitation