PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-02-22 17:38:56 +01:00
Code Releases Activity

fix merge conflict

Browse Source
This commit is contained in:
Antoine Marie 2025-02-21 16:00:57 +01:00
commit 10e8b9c885
84 changed files with 2048 additions and 584 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -6,3 +6,8 @@
__pycache__
.ninja_deps
calcs/
calcs/input
calcs/int
calcs/mol
calcs/water.out

3
PyDuck.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python
#!/usr/bin/env python3
import os, sys
import argparse
@ -163,7 +163,6 @@ if print_2e:
output_file_path = working_dir + '/int/ERI_chem.bin'
subprocess.call(['rm', '-f', output_file_path])
eri_ao = mol.intor('int2e', aosym='s8')
print(eri_ao.shape)
f = open(output_file_path, 'w')
eri_ao.tofile(output_file_path)
f.close()

4
input/methods.default
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF
F F F F
# RHF UHF GHF ROHF HFB
F F F F F
# MP2 MP3
F F
# CCD pCCD DCD CCSD CCSD(T)

4
input/options.default
View File

@ -15,4 +15,6 @@
# ACFDT: AC Kx XBS
F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
F F F F T
F F F F T
# HFB: temperature sigma chem_pot_HF restart_HFB
0.05 1.00 T F

7
mol/CC3_AVTZ/CH2CO.xyz Normal file
View File

@ -0,0 +1,7 @@
5
C 0.00000000 0.00000000 -1.29547953
C 0.00000000 0.00000000 0.01851350
O 0.00000000 0.00000000 1.18357846
H 0.00000000 0.93893013 -1.81881376
H 0.00000000 -0.93893013 -1.81881376

7
mol/CC3_AVTZ/CH2NN.xyz Normal file
View File

@ -0,0 +1,7 @@
5
C 0.00000000 0.00000000 -1.22149978
N 0.00000000 0.00000000 0.07650786
N 0.00000000 0.00000000 1.21670126
H 0.00000000 0.95185857 -1.71597520
H 0.00000000 -0.95185857 -1.71597520

6
mol/CC3_AVTZ/CH2S.xyz Normal file
View File

@ -0,0 +1,6 @@
4
C 0.00000000 0.00000000 -1.10427274
S 0.00000000 0.00000000 0.51463116
H 0.00000000 0.91895736 -1.67756323
H 0.00000000 -0.91895736 -1.677563234

8
mol/CC3_AVTZ/CH3CN.xyz Normal file
View File

@ -0,0 +1,8 @@
7
C 0.00000 0.00000 0.00000
C 1.45672 0.00000 0.00000
H 1.82801 1.02228 0.00000
H 1.82801 -0.51114 0.88532
H 1.82801 -0.51114 -0.88532
N -1.15829 0.00000 0.00000

8
mol/CC3_AVTZ/CH3NC.xyz Normal file
View File

@ -0,0 +1,8 @@
6
C -1.17284 0.00000 0.00000
N 0.00000 0.00000 0.00000
C 1.42106 0.00000 0.00000
H 1.78341 1.02508 0.00000
H 1.78341 -0.51254 0.88775
H 1.78341 -0.51254 -0.88775

8
mol/CC3_AVTZ/CH3NO.xyz Normal file
View File

@ -0,0 +1,8 @@
6
C -0.94419297 0.00000000 -0.56740524
N -0.00286683 0.00000000 0.57183096
O 1.15791903 0.00000000 0.22993880
H -0.40928669 0.00000000 -1.51564611
H -1.57415127 0.88267715 -0.45733920
H -1.57415127 -0.88267715 -0.45733920

7
mol/CC3_AVTZ/HNCO.xyz Normal file
View File

@ -0,0 +1,7 @@
4
C 0.00000 0.00000 0.00000
N 1.21534 0.00000 0.00000
H 1.76899 0.83655 0.00000
O -1.16901 0.00000 0.00000

6
mol/CC3_AVTZ/HOF.xyz Normal file
View File

@ -0,0 +1,6 @@
3
H 0.00000 0.00000 0.00000
O 0.96765 0.00000 0.00000
F 1.16140 1.42521 0.00000

6
mol/CC3_AVTZ/NH2F.xyz Normal file
View File

@ -0,0 +1,6 @@
4
F 0.00000 0.00000 0.00000
N 1.42944 0.00000 0.00000
H 1.62525 0.99892 0.00000
H 1.62525 -0.30787 0.95029

5
mol/CC3_AVTZ/NNO.xyz Normal file
View File

@ -0,0 +1,5 @@
3
N 0.00000 0.00000 0.00000
N 1.13153 0.00000 0.00000
O -1.18937 0.00000 0.00000

108
src/AOtoMO/Hartree_matrix_AO_basis.f90
View File

@ -1,4 +1,4 @@
subroutine Hartree_matrix_AO_basis(nBas,P,G,H)
subroutine Hartree_matrix_AO_basis(nBas,P,ERI,H)
! Compute Hartree matrix in the AO basis
@ -9,7 +9,7 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
@ -25,7 +25,7 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H)
do nu=1,nBas
do la=1,nBas
do mu=1,nBas
H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si)
H(mu,nu) = H(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
end do
end do
end do
@ -45,104 +45,127 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
double precision, intent(in) :: ERI_chem(ERI_size)
double precision, intent(out) :: H(nBas,nBas)
integer :: mu, nu, la, si
integer :: nunu, lala, nula, lasi, numu
integer*8 :: mu, nu, la, si, nBas8
integer*8 :: nunu, lala, nula, lasi, numu
integer*8 :: nunu0, lala0
integer*8 :: nunununu, nunulala, nununula, nunulasi
integer*8 :: lalanunu, lasinunu, numulala, lalanumu
integer*8 :: numunula, numulasi, lasinumu, nununumu
integer*8 :: munu0, munu
integer*8 :: sila0, sila
integer*8 :: munulasi0, munulasi
integer*8 :: nunununu0, numunumu0
do nu = 1, nBas
nunu = (nu * (nu - 1)) / 2 + nu
nunununu = (nunu * (nunu - 1)) / 2 + nunu
nBas8 = int(nBas, kind=8)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (nu, la, si, mu, &
!$OMP nunu0, nunu, nula, lala0, lala, lasi, numu, &
!$OMP nunununu0, nunununu, nununula, numulala, numunula, &
!$OMP nunulala, lalanunu, lalanumu, nunulasi, lasinunu, &
!$OMP numunumu0, nununumu, numulasi, lasinumu) &
!$OMP SHARED (nBas8, H, P, ERI_chem)
!$OMP DO
do nu = 1, nBas8
nunu0 = shiftr(nu * (nu - 1), 1)
nunu = nunu0 + nu
nunununu0 = shiftr(nunu * (nunu - 1), 1)
nunununu = nunununu0 + nunu
H(nu,nu) = P(nu,nu) * ERI_chem(nunununu)
do la = 1, nu-1
do la = 1, nu - 1
lala = (la * (la - 1)) / 2 + la
nunulala = (nunu * (nunu - 1)) / 2 + lala
lala0 = shiftr(la * (la - 1), 1)
lala = lala0 + la
nunulala = nunununu0 + lala
H(nu,nu) = H(nu,nu) + P(la,la) * ERI_chem(nunulala)
nula = (nu * (nu - 1)) / 2 + la
nununula = (nunu * (nunu - 1)) / 2 + nula
nula = nunu0 + la
nununula = nunununu0 + nula
H(nu,nu) = H(nu,nu) + 2.d0 * P(la,nu) * ERI_chem(nununula)
do si = 1, la - 1
lasi = (la * (la - 1)) / 2 + si
nunulasi = (nunu * (nunu - 1)) / 2 + lasi
lasi = lala0 + si
nunulasi = nunununu0 + lasi
H(nu,nu) = H(nu,nu) + 2.d0 * P(si,la) * ERI_chem(nunulasi)
enddo
enddo
do la = nu + 1, nBas
do la = nu + 1, nBas8
lala = (la * (la - 1)) / 2 + la
lalanunu = (lala * (lala - 1)) / 2 + nunu
lala0 = shiftr(la * (la - 1), 1)
lala = lala0 + la
lalanunu = shiftr(lala * (lala - 1), 1) + nunu
H(nu,nu) = H(nu,nu) + P(la,la) * ERI_chem(lalanunu)
do si = 1, la - 1
lasi = (la * (la - 1)) / 2 + si
lasinunu = (lasi * (lasi - 1)) / 2 + nunu
lasi = lala0 + si
lasinunu = shiftr(lasi * (lasi - 1), 1) + nunu
H(nu,nu) = H(nu,nu) + 2.d0 * P(si,la) * ERI_chem(lasinunu)
enddo
enddo
do mu = 1, nu - 1
numu = (nu * (nu - 1)) / 2 + mu
nununumu = (nunu * (nunu - 1)) / 2 + numu
numu = nunu0 + mu
numunumu0 = shiftr(numu * (numu - 1), 1)
nununumu = nunununu0 + numu
H(mu,nu) = p(nu,nu) * ERI_chem(nununumu)
do la = 1, nu - 1
lala = (la * (la - 1)) / 2 + la
numulala = (numu * (numu - 1)) / 2 + lala
lala = shiftr(la * (la - 1), 1) + la
numulala = numunumu0 + lala
H(mu,nu) = H(mu,nu) + p(la,la) * ERI_chem(numulala)
enddo
do la = nu + 1, nBas
lala = (la * (la - 1)) / 2 + la
lalanumu = (lala * (lala - 1)) / 2 + numu
do la = nu + 1, nBas8
lala = shiftr(la * (la - 1), 1) + la
lalanumu = shiftr(lala * (lala - 1), 1) + numu
H(mu,nu) = H(mu,nu) + p(la,la) * ERI_chem(lalanumu)
enddo
do la = 1, mu
nula = (nu * (nu - 1)) / 2 + la
numunula = (numu * (numu - 1)) / 2 + nula
nula = nunu0 + la
numunula = numunumu0 + nula
H(mu,nu) = H(mu,nu) + 2.d0 * P(la,nu) * ERI_chem(numunula)
enddo
do la = mu + 1, nu - 1
nula = (nu * (nu - 1)) / 2 + la
numunula = (nula * (nula - 1)) / 2 + numu
nula = nunu0 + la
numunula = shiftr(nula * (nula - 1), 1) + numu
H(mu,nu) = H(mu,nu) + 2.d0 * P(la,nu) * ERI_chem(numunula)
enddo
do la = 2, nu - 1
lala0 = shiftr(la * (la - 1), 1)
do si = 1, la - 1
lasi = (la * (la - 1)) / 2 + si
numulasi = (numu * (numu - 1)) / 2 + lasi
lasi = lala0 + si
numulasi = numunumu0 + lasi
H(mu,nu) = H(mu,nu) + 2.d0 * P(si,la) * ERI_chem(numulasi)
enddo
enddo
do la = nu + 1, nBas
do la = nu + 1, nBas8
lala0 = shiftr(la * (la - 1), 1)
do si = 1, la - 1
lasi = (la * (la - 1)) / 2 + si
lasinumu = (lasi * (lasi - 1)) / 2 + numu
lasi = lala0 + si
lasinumu = shiftr(lasi * (lasi - 1), 1) + numu
H(mu,nu) = H(mu,nu) + 2.d0 * P(si,la) * ERI_chem(lasinumu)
enddo
enddo
enddo ! mu
enddo ! nu
!$OMP END DO
!$OMP END PARALLEL
do nu = 1, nBas
do mu = nu+1, nBas
do nu = 1, nBas8
do mu = nu+1, nBas8
H(mu,nu) = H(nu,mu)
enddo
enddo
@ -152,3 +175,4 @@ end subroutine
! ---

38
src/AOtoMO/anomalous_matrix_AO_basis.f90 Normal file
View File

@ -0,0 +1,38 @@
subroutine anomalous_matrix_AO_basis(nBas,sigma,Pa,ERI,L)
! Compute anomalous L matrix in the AO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: sigma
double precision,intent(in) :: Pa(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: L(nBas,nBas)
L(:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do la=1,nBas
do mu=1,nBas
L(mu,nu) = L(mu,nu) + sigma*Pa(la,si)*ERI(la,si,mu,nu)
end do
end do
end do
end do
end subroutine
! ---

201
src/AOtoMO/exchange_matrix_AO_basis.f90
View File

@ -19,7 +19,8 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K)
double precision,intent(out) :: K(nBas,nBas)
K = 0d0
K(:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do la=1,nBas
@ -44,55 +45,68 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
double precision, intent(in) :: ERI_chem(ERI_size)
double precision, intent(out) :: K(nBas,nBas)
integer :: mu, nu, la, si
integer :: nunu, nula, lanu, lasi, nusi, sinu
integer :: numu, mumu, mula, lamu, musi, simu
integer*8 :: mu, nu, la, si, nBas8
integer*8 :: nunu, nula, lanu, lasi, nusi, sinu
integer*8 :: numu, mumu, mula, lamu, musi, simu
integer*8 :: nunu0, lala0, mumu0
integer*8 :: nunununu, nulanula, lanulanu, nulanusi
integer*8 :: munulasi, lanunusi, lanusinu, numumumu
integer*8 :: nulamula, nulalamu, lanulamu, nulamusi
integer*8 :: nulasimu, lanumusi, lanusimu
integer*8 :: nulasimu, lanumusi, lanusimu, simunula
integer*8 :: simulanu, nulanula0, lanulanu0
integer*8, external :: Yoshimine_4ind
nBas8 = int(nBas, kind=8)
do nu = 1, nBas
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (nu, si, la, mu, &
!$OMP nunu0, nunu, lanu, numu, mumu0, mumu, simu, lala0, nula, &
!$OMP nunununu, nulanula, lanulanu, lanulanu0, nulanula0, &
!$OMP nulanusi, lanulamu, lanunusi, lanusinu , numumumu, &
!$OMP nulamula, nulalamu, lanumusi, lanusimu, nulamusi, &
!$OMP nulasimu, simunula, simulanu) &
!$OMP SHARED (nBas8, P, ERI_chem, K)
!$OMP DO
do nu = 1, nBas8
nunu = (nu * (nu - 1)) / 2 + nu
nunununu = (nunu * (nunu - 1)) / 2 + nunu
nunu0 = shiftr(nu * (nu - 1), 1)
nunu = nunu0 + nu
nunununu = shiftr(nunu * (nunu - 1), 1) + nunu
K(nu,nu) = -P(nu,nu) * ERI_chem(nunununu)
do la = 1, nu - 1
nula = (nu * (nu - 1)) / 2 + la
nulanula = (nula * (nula - 1)) / 2 + nula
nula = nunu0 + la
nulanula = shiftr(nula * (nula - 1), 1) + nula
K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(nulanula)
enddo
do la = nu + 1, nBas
lanu = (la * (la - 1)) / 2 + nu
lanulanu = (lanu * (lanu - 1)) / 2 + lanu
do la = nu + 1, nBas8
lanu = shiftr(la * (la - 1), 1) + nu
lanulanu = shiftr(lanu * (lanu - 1), 1) + lanu
K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(lanulanu)
enddo
do la = 1, nu
nula = (nu * (nu - 1)) / 2 + la
nula = nunu0 + la
nulanula0 = shiftr(nula * (nula - 1), 1)
do si = 1, la - 1
nusi = (nu * (nu - 1)) / 2 + si
nulanusi = (nula * (nula - 1)) / 2 + nusi
nulanusi = nulanula0 + nunu0 + si
K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(nulanusi)
enddo
enddo
do la = nu + 1, nBas
lanu = (la * (la - 1)) / 2 + nu
do la = nu + 1, nBas8
lanu = shiftr(la * (la - 1), 1) + nu
lanulanu0 = shiftr(lanu * (lanu - 1), 1)
do si = 1, nu
nusi = (nu * (nu - 1)) / 2 + si
lanunusi = (lanu * (lanu - 1)) / 2 + nusi
lanunusi = lanulanu0 + nunu0 + si
K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanunusi)
enddo
do si = nu + 1, la - 1
sinu = (si * (si - 1)) / 2 + nu
lanusinu = (lanu * (lanu - 1)) / 2 + sinu
lanusinu = lanulanu0 + shiftr(si * (si - 1), 1) + nu
K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanusinu)
enddo
enddo
@ -100,115 +114,109 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
do mu = 1, nu - 1
numu = (nu * (nu - 1)) / 2 + mu
mumu = (mu * (mu - 1)) / 2 + mu
numumumu = (numu * (numu - 1)) / 2 + mumu
numu = nunu0 + mu
mumu0 = shiftr(mu * (mu - 1), 1)
mumu = mumu0 + mu
numumumu = shiftr(numu * (numu - 1), 1) + mumu
K(mu,nu) = - P(mu,mu) * ERI_chem(numumumu)
do la = 1, mu - 1
mula = (mu * (mu - 1)) / 2 + la
nula = (nu * (nu - 1)) / 2 + la
nulamula = (nula * (nula - 1)) / 2 + mula
nula = nunu0 + la
nulamula = shiftr(nula * (nula - 1), 1) + mumu0 + la
K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulamula)
enddo
do la = mu + 1, nu
lamu = (la * (la - 1)) / 2 + mu
nula = (nu * (nu - 1)) / 2 + la
nulalamu = (nula * (nula - 1)) / 2 + lamu
nula = nunu0 + la
nulalamu = shiftr(nula * (nula - 1), 1) + shiftr(la * (la - 1), 1) + mu
K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulalamu)
enddo
do la = nu + 1, nBas
lamu = (la * (la - 1)) / 2 + mu
lanu = (la * (la - 1)) / 2 + nu
lanulamu = (lanu * (lanu - 1)) / 2 + lamu
do la = nu + 1, nBas8
lala0 = shiftr(la * (la - 1), 1)
lanu = lala0 + nu
lanulamu = shiftr(lanu * (lanu - 1), 1) + lala0 + mu
K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(lanulamu)
enddo
do la = 1, mu
nula = (nu * (nu - 1)) / 2 + la
nula = nunu0 + la
nulanula0 = shiftr(nula * (nula - 1), 1)
do si = 1, la - 1
musi = (mu * (mu - 1)) / 2 + si
nulamusi = (nula * (nula - 1)) / 2 + musi
nulamusi = nulanula0 + mumu0 + si
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi)
enddo
enddo
do la = mu + 1, nu
nula = (nu * (nu - 1)) / 2 + la
nula = nunu0 + la
nulanula0 = shiftr(nula * (nula - 1), 1)
do si = 1, mu
musi = (mu * (mu - 1)) / 2 + si
nulamusi = (nula * (nula - 1)) / 2 + musi
nulamusi = nulanula0 + mumu0 + si
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi)
enddo
do si = mu + 1, la - 1
simu = (si * (si - 1)) / 2 + mu
nulasimu = (nula * (nula - 1)) / 2 + simu
nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu)
enddo
enddo
do la = nu + 1, nBas
lanu = (la * (la - 1)) / 2 + nu
do la = nu + 1, nBas8
lanu = shiftr(la * (la - 1), 1) + nu
lanulanu0 = shiftr(lanu * (lanu - 1), 1)
do si = 1, mu
musi = (mu * (mu - 1)) / 2 + si
lanumusi = (lanu * (lanu - 1)) / 2 + musi
lanumusi = lanulanu0 + mumu0 + si
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanumusi)
enddo
do si = mu + 1, la - 1
simu = (si * (si - 1)) / 2 + mu
lanusimu = (lanu * (lanu - 1)) / 2 + simu
lanusimu = lanulanu0 + shiftr(si * (si - 1), 1) + mu
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanusimu)
enddo
enddo
!TODO
! do la = 1, mu
! nula = (nu * (nu - 1)) / 2 + la
! do si = la + 1, mu
! musi = (mu * (mu - 1)) / 2 + si
! nulamusi = (nula * (nula - 1)) / 2 + musi
! !nulamusi = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi)
! enddo
! do si = mu + 1, nBas
! simu = (si * (si - 1)) / 2 + mu
! nulasimu = (nula * (nula - 1)) / 2 + simu
! !nulasimu = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu)
! enddo
! enddo
! do la = mu + 1, nu
! nula = (nu * (nu - 1)) / 2 + la
! do si = la + 1, nu
! simu = (si * (si - 1)) / 2 + mu
! nulasimu = (nula * (nula - 1)) / 2 + simu
! !nulasimu = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu)
! enddo
! do si = nu + 1, nBas
! simu = (si * (si - 1)) / 2 + mu
! munulasi = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
! enddo
! enddo
! do la = nu + 1, nBas
! lanu = (la * (la - 1)) / 2 + nu
! do si = la + 1, mu
! simu = (si * (si - 1)) / 2 + mu
! munulasi = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
! enddo
! do si = mu + 1, nBas
! musi = (mu * (mu - 1)) / 2 + si
! munulasi = Yoshimine_4ind(nu, la, si, mu)
! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
! enddo
! enddo
do la = 1, mu
nula = nunu0 + la
nulanula0 = shiftr(nula * (nula - 1), 1)
do si = la + 1, mu
nulamusi = nulanula0 + mumu0 + si
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi)
enddo
do si = mu + 1, nu - 1
nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu)
enddo
do si = nu, nBas8
simu = shiftr(si * (si - 1), 1) + mu
simunula = shiftr(simu * (simu - 1), 1) + nula
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simunula)
enddo
enddo
do la = mu + 1, nu
nula = nunu0 + la
nulanula0 = shiftr(nula * (nula - 1), 1)
do si = la + 1, nu
nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu)
enddo
do si = nu + 1, nBas8
simu = shiftr(si * (si - 1), 1) + mu
simunula = shiftr(simu * (simu - 1), 1) + nula
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simunula)
enddo
enddo
do la = nu + 1, nBas8
lanu = shiftr(la * (la - 1), 1) + nu
do si = la + 1, nBas8
simu = shiftr(si * (si - 1), 1) + mu
simulanu = shiftr(simu * (simu - 1), 1) + lanu
K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simulanu)
enddo
enddo
enddo ! mu
enddo ! nu
!$OMP END DO
!$OMP END PARALLEL
do nu = 1, nBas
do mu = nu+1, nBas
do nu = 1, nBas8
do mu = nu+1, nBas8
K(mu,nu) = K(nu,mu)
enddo
enddo
@ -218,3 +226,6 @@ end subroutine
! ---

144
src/CC/EE_EOM_CCD_1h1p.f90
View File

@ -28,10 +28,22 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
double precision,allocatable :: Wovvo(:,:,:,:)
double precision,allocatable :: H(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: VL(:,:)
double precision,allocatable :: VR(:,:)
double precision,allocatable :: Leom(:,:,:)
double precision,allocatable :: Reom(:,:,:)
integer :: nstate,m
double precision :: Ex,tmp
integer,allocatable :: order(:)
double precision,allocatable :: rdm1_oo(:,:)
double precision,allocatable :: rdm1_vv(:,:)
double precision,allocatable :: rdm2_oovv(:,:,:,:)
double precision,allocatable :: rdm2_ovvo(:,:,:,:)
! Hello world
write(*,*)
@ -46,10 +58,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
! Memory allocation
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),Z(nS,nS))
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS))
allocate(order(nS))
! Form one-body terms
do a=1,nV-nR
@ -132,19 +143,142 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
! Diagonalize EOM Hamiltonian
allocate(VL(nS,nS),VR(nS,nS))
if(nS > 0) then
call diagonalize_general_matrix(nS,H,Om,Z)
call diagonalize_general_matrix_LR(nS,H,Om,VL,VR)
do ia=1,nS
order(ia) = ia
end do
call quick_sort(Om,order,nS)
call set_order(Z,order,nS,nS)
call set_order_LR(VL,VR,order,nS,nS)
call print_excitation_energies('EE-EOM-CCD','spinorbital',nS,Om)
! write(*,*) 'Right Eigenvectors'
! call matout(nS,nS,VR)
! call matout(nS,3,VR(:,1:3))
end if
allocate(Leom(nO,nV,nS),Reom(nO,nV,nS))
do m=1,nS
ia = 0
do i=1,nO
do a=1,nV
ia = ia + 1
Leom(i,a,m) = VL(ia,m)
Reom(i,a,m) = VR(ia,m)
end do
end do
end do
deallocate(VL,VR)
!------------------------------------------------------------------------
! EOM section
!------------------------------------------------------------------------
allocate(rdm1_oo(nO,nO),rdm1_vv(nV,nV))
allocate(rdm2_oovv(nO,nO,nV,nV),rdm2_ovvo(nO,nV,nV,nO))
nstate = 1
tmp = 0d0
do i=1,nO
do a=1,nV
tmp = tmp + Leom(i,a,nstate)*Reom(i,a,nstate)
end do
end do
print*,tmp
rdm1_oo(:,:) = 0d0
do i=1,nO
do j=1,nO
do c=1,nV
rdm1_oo(i,j) = rdm1_oo(i,j) - Reom(i,c,nstate)*Leom(j,c,nstate)
end do
end do
end do
rdm1_vv(:,:) = 0d0
do a=1,nV
do b=1,nV
do k=1,nO
rdm1_vv(a,b) = rdm1_vv(a,b) + Reom(k,b,nstate)*Leom(k,a,nstate)
end do
end do
end do
rdm2_ovvo(:,:,:,:) = 0d0
do i=1,nO
do a=1,nV
do b=1,nV
do j=1,nO
rdm2_ovvo(i,a,b,j) = Reom(i,b,nstate)*Leom(j,a,nstate)
end do
end do
end do
end do
rdm2_oovv(:,:,:,:) = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do k=1,nO
do c=1,nV
rdm2_oovv(i,j,a,b) = rdm2_oovv(i,j,a,b) &
+ Reom(j,b,nstate)*t(k,i,c,a)*Leom(k,c,nstate) &
- Reom(i,b,nstate)*t(k,j,c,a)*Leom(k,c,nstate) &
- Reom(j,a,nstate)*t(k,i,c,b)*Leom(k,c,nstate) &
+ Reom(i,a,nstate)*t(k,j,c,b)*Leom(k,c,nstate)
end do
end do
end do
end do
end do
end do
Ex = 0d0
do i=1,nO
Ex = Ex + rdm1_oo(i,i)*eO(i)
end do
do a=1,nV
Ex = Ex + rdm1_vv(a,a)*eV(a)
end do
do i=1,nO
do a=1,nV
do b=1,nV
do j=1,nO
Ex = Ex + rdm2_ovvo(i,a,b,j)*OVVO(i,a,b,j) + 0.25d0*rdm2_oovv(i,j,a,b)*OOVV(i,j,a,b)
end do
end do
end do
end do
print*,'Ex = ',Ex
print*,'Om = ',Om(nstate)
end subroutine

1
src/CC/rCCD.f90
View File

@ -222,7 +222,6 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
if(dotest) then
call dump_test_value('R','rCCD correlation energy',EcCC)

4
src/CI/RCIS.f90
View File

@ -48,7 +48,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
if(singlet) then
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (singlet state)'
@ -84,7 +84,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
if(triplet) then
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (triplet state)'

2
src/GF/GGF2_phBSE.f90
View File

@ -57,7 +57,7 @@ subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,E
! Compute phBSE@GF2 excitation energies
call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
call phGLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
call print_excitation_energies('phBSE@GGF2','generalized',nS,OmBSE)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)

2
src/GF/GGF2_ppBSE.f90
View File

@ -72,7 +72,7 @@ subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBS
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)

12
src/GF/RGF2_phBSE.f90
View File

@ -54,8 +54,8 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
ispin = 1
EcBSE(ispin) = 0d0
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
! Compute static kernel
@ -67,7 +67,7 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
! Compute phBSE@GF2 excitation energies
call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
call phRLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
call print_excitation_energies('phBSE@GF2','singlet',nS,OmBSE)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
@ -87,8 +87,8 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
ispin = 2
EcBSE(ispin) = 0d0
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
! Compute static kernel
@ -100,7 +100,7 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
! Compute phBSE@GF2 excitation energies
call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
call phRLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
call print_excitation_energies('phBSE@GF2','triplet',nS,OmBSE)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)

16
src/GF/RGF2_ppBSE.f90
View File

@ -78,15 +78,15 @@ subroutine RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dip
call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
@ -131,15 +131,15 @@ subroutine RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dip
call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)

12
src/GT/RG0T0eh.f90
View File

@ -109,10 +109,10 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
ispin = 2
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
@ -158,10 +158,10 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------!
! Dump results !

32
src/GT/RG0T0pp.f90
View File

@ -115,11 +115,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppRLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -136,11 +136,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppRLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -207,11 +207,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppRLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
deallocate(Bpp,Cpp,Dpp)
@ -220,11 +220,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppRLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
deallocate(Bpp,Cpp,Dpp)

44
src/GT/RGTpp_phACFTD.f90
View File

@ -140,11 +140,11 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppRLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp)
call ppRLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -158,11 +158,11 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppRLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp)
call ppRLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -173,13 +173,13 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
end if
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
Aph(:,:) = Aph(:,:) - TAt(:,:) - TAs(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) - TBt(:,:) - TBs(:,:)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
@ -224,11 +224,11 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppRLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp)
call ppRLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -242,11 +242,11 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppRLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp)
call ppRLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
@ -257,13 +257,13 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
end if
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
Aph(:,:) = Aph(:,:) + TAs(:,:) - TAt(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBs(:,:) - TBt(:,:)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))

28
src/GT/RGTpp_phBSE.f90
View File

@ -97,11 +97,11 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
@ -116,11 +116,11 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
@ -137,13 +137,13 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
! Compute BSE singlet excitation energies
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
Aph(:,:) = Aph(:,:) + TAt(:,:) ! TAt(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) ! TBt(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GTpp','singlet',nS,OmBSE)
call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
@ -167,13 +167,13 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
! Compute BSE triplet excitation energies
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
Aph(:,:) = Aph(:,:) + 1d0*TAt(:,:) - TAs(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + 1d0*TBt(:,:) - TBs(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GTpp','triplet',nS,OmBSE)
call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)

32
src/GT/RGTpp_ppBSE.f90
View File

@ -65,14 +65,14 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppRLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs))
allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs))
@ -87,14 +87,14 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppRLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt))
allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
allocate(rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
@ -146,9 +146,9 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
allocate(KB_sta(nVVs,nOOs),KC_sta(nVVs,nVVs),KD_sta(nOOs,nOOs))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KB_sta)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KC_sta)
@ -158,7 +158,7 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
call ppRLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
@ -184,9 +184,9 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
allocate(KB_sta(nVVt,nOOt),KC_sta(nVVt,nVVt),KD_sta(nOOt,nOOt))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,Taaaa,Tabab,Tbaab,KB_sta)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,Taaaa,Tabab,Tbaab,KC_sta)
@ -196,7 +196,7 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin))
call ppRLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin))
call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)

6
src/GT/evRGTeh.f90
View File

@ -123,10 +123,10 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute screening
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phrLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)

16
src/GT/evRGTpp.f90
View File

@ -135,11 +135,11 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
@ -153,11 +153,11 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)

25
src/GT/qsRGTeh.f90
View File

@ -197,10 +197,10 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! Compute linear response
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
@ -221,10 +221,6 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:)
if(nBas .ne. nOrb) then
call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
endif
! Compute commutator and convergence criteria
@ -241,17 +237,10 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! Diagonalize Hamiltonian in AO basis
if(nBas .eq. nOrb) then
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eGT)
c = matmul(X, cp)
else
Fp = matmul(transpose(c), matmul(F, c))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eGT)
c = matmul(c, cp)
endif
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eGT)
c = matmul(X, cp)
! Compute new density matrix in the AO basis

16
src/GT/qsRGTpp.f90
View File

@ -198,11 +198,11 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
@ -210,11 +210,11 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)

16
src/GT/ufRG0T0pp.f90
View File

@ -115,11 +115,11 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppRLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s(:))
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s(:))
@ -138,11 +138,11 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppRLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOt,Om2t)

12
src/GW/RG0W0.f90
View File

@ -103,10 +103,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Compute screening !
!-------------------!
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
@ -160,10 +160,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Compute the RPA correlation energy
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------!
! Dump results !

30
src/GW/RGW_phACFDT.f90
View File

@ -93,10 +93,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
isp_W = 1
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
@ -123,10 +123,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
if(doXBS) then
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
@ -134,13 +134,13 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
end if
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
Aph(:,:) = Aph(:,:) + KA(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
@ -180,10 +180,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
if(doXBS) then
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
@ -191,13 +191,13 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
end if
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
Aph(:,:) = Aph(:,:) + KA(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))

18
src/GW/RGW_phBSE.f90
View File

@ -83,10 +83,10 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
isp_W = 1
EcRPA = 0d0
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
@ -102,8 +102,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
! Compute BSE excitation energies
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
! Second-order BSE static kernel
@ -127,7 +127,7 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
Aph(:,:) = Aph(:,:) + KA_sta(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GW@RHF','singlet',nS,OmBSE)
call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
@ -160,13 +160,13 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
! Compute BSE excitation energies
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
Aph(:,:) = Aph(:,:) + KA_sta(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GW@RHF','triplet',nS,OmBSE)
call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)

25
src/GW/RGW_ppBSE.f90
View File

@ -98,11 +98,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
Aph(nS,nS),Bph(nS,nS))
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
! call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
@ -143,17 +142,15 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
endif
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
if(.not.TDA) then
call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
endif
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
deallocate(Bpp,Cpp,Dpp)
!
@ -265,15 +262,15 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
if(.not.TDA) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
call ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
deallocate(Bpp,Cpp,Dpp)
! ---

6
src/GW/ccRG0W0.f90
View File

@ -93,10 +93,10 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI,ENuc,ERHF
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phRLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)

6
src/GW/evRGW.f90
View File

@ -120,10 +120,10 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Compute screening
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute spectral weights

6
src/GW/qsRGW.f90
View File

@ -202,10 +202,10 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Compute linear response
call phLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phRLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phRLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)

6
src/GW/ufRG0W0.f90
View File

@ -96,10 +96,10 @@ subroutine ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phRLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)

6
src/GW/ufRGW.f90
View File

@ -96,10 +96,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ispin = 1
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)

373
src/HF/HFB.f90 Normal file
View File

@ -0,0 +1,373 @@
subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,Panom,F,Delta, &
temperature,sigma,chem_pot_hf,restart_hfb)
! Perform Hartree-Fock Bogoliubov calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature,sigma
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
logical :: chem_pot_hf
logical :: restart_hfb
integer :: nBas2
integer :: iorb
integer :: nSCF
integer :: nOrb2
integer :: nBas2_Sq
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: EL
double precision :: chem_pot
double precision :: dipole(ncart)
double precision :: Conv
double precision :: rcond
double precision :: trace_1rdm
double precision :: thrs_N
double precision :: norm_anom
double precision,external :: trace_matrix
double precision,allocatable :: eigVAL(:)
double precision,allocatable :: Occ(:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: H_HFB_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: eigVEC(:,:)
double precision,allocatable :: H_HFB(:,:)
double precision,allocatable :: R(:,:)
double precision,allocatable :: err_ao(:,:)
double precision,allocatable :: S_ao(:,:)
double precision,allocatable :: X_ao(:,:)
double precision,allocatable :: R_ao_old(:,:)
double precision,allocatable :: H_HFB_ao(:,:)
! Output variables
double precision,intent(out) :: EHFB
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: Panom(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
double precision,intent(out) :: Delta(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'*****************************'
write(*,*)'* HF Bogoliubov Calculation *'
write(*,*)'*****************************'
write(*,*)
! Useful quantities
nOrb2 = nOrb+nOrb
nBas2 = nBas+nBas
nBas2_Sq = nBas2*nBas2
! Memory allocation
allocate(Occ(nOrb))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(eigVEC(nOrb2,nOrb2))
allocate(H_HFB(nOrb2,nOrb2))
allocate(R(nOrb2,nOrb2))
allocate(eigVAL(nOrb2))
allocate(err_ao(nBas2,nBas2))
allocate(S_ao(nBas2,nBas2))
allocate(X_ao(nBas2,nOrb2))
allocate(R_ao_old(nBas2,nBas2))
allocate(H_HFB_ao(nBas2,nBas2))
allocate(err_diis(nBas2_Sq,max_diis))
allocate(H_HFB_diis(nBas2_Sq,max_diis))
! Guess chem. pot.
chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
! Initialization
thrs_N = 1d-8
n_diis = 0
H_HFB_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
P(:,:) = matmul(c(:,1:nO), transpose(c(:,1:nO)))
Panom(:,:) = 0d0
! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot
if(abs(temperature)>1d-4) then
Occ(:) = 0d0
Occ(1:nO) = 1d0
call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
if(chem_pot_hf) chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
P(:,:) = 0d0
Panom(:,:) = 0d0
do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
endif
! Read restart file
if(restart_hfb) then
call read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
P(:,:) = 0d0
Panom(:,:) = 0d0
do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
endif
P(:,:) = 2d0*P(:,:)
S_ao(:,:) = 0d0
S_ao(1:nBas ,1:nBas ) = S(1:nBas,1:nBas)
S_ao(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas)
X_ao(:,:) = 0d0
X_ao(1:nBas ,1:nOrb ) = X(1:nBas,1:nOrb)
X_ao(nBas+1:nBas2,nOrb+1:nOrb2) = X(1:nBas,1:nOrb)
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*) 'Enterning HFB SCF procedure'
write(*,*)
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock and Delta matrices
call Hartree_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
call anomalous_matrix_AO_basis(nBas,sigma,Panom,ERI,Delta)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
! Diagonalize H_HFB matrix
H_HFB(:,:) = 0d0
H_HFB(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X))
H_HFB(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_HFB(1:nOrb,1:nOrb)
H_HFB(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X))
H_HFB(nOrb+1:nOrb2,1:nOrb ) = H_HFB(1:nOrb,nOrb+1:nOrb2)
eigVEC(:,:) = H_HFB(:,:)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
! Build R
R(:,:) = 0d0
do iorb=1,nOrb
R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb)))
enddo
trace_1rdm = 0d0
do iorb=1,nOrb
trace_1rdm = trace_1rdm+R(iorb,iorb)
enddo
! Adjust the chemical potential
if( abs(trace_1rdm-nO) > thrs_N ) &
call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL)
! DIIS extrapolation
if(max_diis > 1 .and. nSCF>1) then
write(*,*) ' Doing DIIS'
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
H_HFB_ao(:,:) = 0d0
H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas)
H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas)
err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao)
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas2_Sq,nBas2_Sq,n_diis,err_diis,H_HFB_diis,err_ao,H_HFB_ao)
H_HFB = matmul(transpose(X_ao),matmul(H_HFB_ao,X_ao))
eigVEC(:,:) = H_HFB(:,:)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
! Build R and check trace
trace_1rdm = 0d0
R(:,:) = 0d0
do iorb=1,nOrb
R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb)))
enddo
do iorb=1,nOrb
trace_1rdm = trace_1rdm + R(iorb,iorb)
enddo
! Adjust the chemical potential
if( abs(trace_1rdm-nO) > thrs_N ) &
call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL)
end if
! Extract P and Panom from R
P(:,:) = 0d0
Panom(:,:) = 0d0
P(:,:) = 2d0*matmul(X,matmul(R(1:nOrb,1:nOrb),transpose(X)))
Panom(:,:) = matmul(X,matmul(R(1:nOrb,nOrb+1:nOrb2),transpose(X)))
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Anomalous energy
EL = trace_matrix(nBas,matmul(Panom,Delta))
! Total energy
EHFB = ET + EV + EJ + EK + EL
! Check convergence
if(nSCF > 1) then
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
H_HFB_ao(:,:) = 0d0
H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas)
H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas)
err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao)
Conv = maxval(abs(err_ao))
endif
! Update R_old
R_ao_old(:,:) = 0d0
R_ao_old(1:nBas ,1:nBas ) = 0.5d0*P(1:nBas,1:nBas)
R_ao_old(nBas+1:nBas2,nBas+1:nBas2) = matmul(X(1:nBas,1:nOrb), transpose(X(1:nBas,1:nOrb)))-0.5d0*P(1:nBas,1:nBas)
R_ao_old(1:nBas ,nBas+1:nBas2) = Panom(1:nBas,1:nBas)
R_ao_old(nBas+1:nBas2,1:nBas ) = Panom(1:nBas,1:nBas)
! Dump results
write(*,*)'-----------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') &
'|','#','|','E(HFB)','|','EJ(HFB)','|','EK(HFB)','|','EL(HFB)','|','Conv','|'
write(*,*)'-----------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1XF16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',EHFB + ENuc,'|',EJ,'|',EK,'|',EL,'|',Conv,'|'
write(*,*)'-----------------------------------------------------------------------------------------------'
write(*,*)
end do
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
stop
end if
! Compute dipole moments, occupation numbers and || Anomalous density||
! also print the restart file
deallocate(eigVEC,eigVAL)
allocate(eigVEC(nOrb,nOrb),eigVAL(nOrb))
eigVEC(1:nOrb,1:nOrb) = 0d0
eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb)
call diagonalize_matrix(nOrb,eigVEC,eigVAL)
Occ(1:nOrb) = eigVAL(1:nOrb)
c = matmul(X,eigVEC)
norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2)))
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot)
call print_HFB(nBas,nOrb,nO,norm_anom,Occ,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole)
! Testing zone
if(dotest) then
call dump_test_value('R','HFB energy',EHFB)
call dump_test_value('R','HFB dipole moment',norm2(dipole))
end if
! Memory deallocation
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
end subroutine

52
src/HF/MOM_overlap.f90
View File

@ -1,52 +0,0 @@
subroutine MOM_overlap(nBas,nO,S,cG,c,ON)
! Compute overlap between old and new MO coefficients
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: S(nBas,nBas),cG(nBas,nBas),c(nBas,nBas)
! Local variables
integer :: i,j,ploc
double precision,allocatable :: Ov(:,:),pOv(:)
! Output variables
double precision,intent(inout):: ON(nBas)
allocate(Ov(nBas,nBas),pOv(nBas))
Ov = matmul(transpose(cG),matmul(S,c))
pOv(:) = 0d0
do i=1,nBas
do j=1,nBas
pOv(j) = pOv(j) + Ov(i,j)**2
end do
end do
pOv(:) = sqrt(pOV(:))
! print*,'--- MOM overlap ---'
! call matout(nBas,1,pOv)
ON(:) = 0d0
do i=1,nO
ploc = maxloc(pOv,nBas)
ON(ploc) = 1d0
pOv(ploc) = 0d0
end do
! print*,'--- Occupation numbers ---'
! call matout(nBas,1,ON)
end subroutine

15
src/HF/RHF.f90
View File

@ -127,10 +127,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! if(nBas .ne. nOrb) then
! call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
! call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
! endif
! Check convergence
@ -174,21 +170,14 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Diagonalize Fock matrix
! if(nBas .eq. nOrb) then
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
! else
! Fp = matmul(transpose(c),matmul(F,c))
! cp(:,:) = Fp(:,:)
! call diagonalize_matrix(nOrb,cp,eHF)
! c = matmul(c,cp)
! endif
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
! c(1,1), nBas, c(1,1), nBas, &
! 0.d0, P(1,1), nBas)
@ -236,6 +225,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
end if
! Memory deallocation
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
end subroutine

152
src/HF/RHF_hpc.f90
View File

@ -39,6 +39,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
integer :: nBas_Sq
integer :: n_diis
integer*8 :: ERI_size
double precision :: t1, t2
double precision :: diff, diff_loc
double precision :: ET
double precision :: EV
@ -57,6 +58,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ERI_chem(:)
double precision,allocatable :: ERI_phys(:,:,:,:), J_deb(:,:), K_deb(:,:)
double precision,allocatable :: tmp1(:,:), FX(:,:)
! Output variables
@ -92,6 +94,9 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
allocate(tmp1(nBas,nBas))
allocate(FX(nBas,nOrb))
! Guess coefficients and density matrix
call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
@ -99,52 +104,60 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)
ERI_size = (nBas * (nBas + 1)) / 2
ERI_size = (ERI_size * (ERI_size + 1)) / 2
ERI_size = shiftr(nBas * (nBas + 1), 1)
ERI_size = shiftr(ERI_size * (ERI_size + 1), 1)
allocate(ERI_chem(ERI_size))
call read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem)
call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J)
allocate(ERI_phys(nBas,nBas,nBas,nBas))
allocate(J_deb(nBas,nBas))
allocate(K_deb(nBas,nBas))
call read_2e_integrals(working_dir, nBas, ERI_phys)
call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
print*, "max error on J = ", maxval(dabs(J - J_deb))
diff = 0.d0
do ii = 1, nBas
do jj = 1, nBas
diff_loc = dabs(J(jj,ii) - J_deb(jj,ii))
if(diff_loc .gt. 1d-12) then
print*, 'error in J on: ', jj, ii
print*, J(jj,ii), J_deb(jj,ii)
stop
endif
diff = diff + diff_loc
enddo
enddo
print*, 'total diff on J = ', diff
call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb)
print*, "max error on K = ", maxval(dabs(K - K_deb))
diff = 0.d0
do ii = 1, nBas
do jj = 1, nBas
diff_loc = dabs(K(jj,ii) - K_deb(jj,ii))
if(diff_loc .gt. 1d-12) then
print*, 'error in K on: ', jj, ii
print*, K(jj,ii), K_deb(jj,ii)
stop
endif
diff = diff + diff_loc
enddo
enddo
print*, 'total diff on K = ', diff
stop
!call wall_time(t1)
!call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J)
!call wall_time(t2)
!print*, " J built in (sec):", (t2-t1)
!call wall_time(t1)
!call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
!call wall_time(t2)
!print*, " K built in (sec):", (t2-t1)
!allocate(ERI_phys(nBas,nBas,nBas,nBas))
!allocate(J_deb(nBas,nBas))
!allocate(K_deb(nBas,nBas))
!call read_2e_integrals(working_dir, nBas, ERI_phys)
!call wall_time(t1)
!call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
!call wall_time(t2)
!print*, " J_deb built in (sec):", (t2-t1)
!call wall_time(t1)
!call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb)
!call wall_time(t2)
!print*, " K_deb built in (sec):", (t2-t1)
!print*, "max error on J = ", maxval(dabs(J - J_deb))
!diff = 0.d0
!do ii = 1, nBas
! do jj = 1, nBas
! diff_loc = dabs(J(jj,ii) - J_deb(jj,ii))
! if(diff_loc .gt. 1d-10) then
! print*, 'error in J on: ', jj, ii
! print*, J(jj,ii), J_deb(jj,ii)
! stop
! endif
! diff = diff + diff_loc
! enddo
!enddo
!print*, 'total diff on J = ', diff
!print*, "max error on K = ", maxval(dabs(K - K_deb))
!diff = 0.d0
!do ii = 1, nBas
! do jj = 1, nBas
! diff_loc = dabs(K(jj,ii) - K_deb(jj,ii))
! if(diff_loc .gt. 1d-10) then
! print*, 'error in K on: ', jj, ii
! print*, K(jj,ii), K_deb(jj,ii)
! stop
! endif
! diff = diff + diff_loc
! enddo
!enddo
!print*, 'total diff on K = ', diff
!stop
! Initialization
@ -173,25 +186,46 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
nSCF = nSCF + 1
! Build Fock matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_phys,J)
call exchange_matrix_AO_basis(nBas,P,ERI_phys,K)
call Hartree_matrix_AO_basis_hpc (nBas, ERI_size, P(1,1), ERI_chem(1), J(1,1))
call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P(1,1), ERI_chem(1), K(1,1))
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! Check convergence
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
S(1,1), nBas, P(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
tmp1(1,1), nBas, F(1,1), nBas, &
0.d0, err(1,1), nBas)
call dgemm("N", "T", nBas, nBas, nBas, 1.d0, &
F(1,1), nBas, tmp1(1,1), nBas, &
-1.d0, err(1,1), nBas)
!err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy
ET = trace_matrix(nBas, matmul(P, T))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, T(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
ET = trace_matrix(nBas, tmp1(1,1))
! Potential energy
EV = trace_matrix(nBas, matmul(P, V))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, V(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EV = trace_matrix(nBas, tmp1(1,1))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, J(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EJ = 0.5d0*trace_matrix(nBas, tmp1(1,1))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas, matmul(P, K))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, K(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EK = 0.25d0*trace_matrix(nBas, tmp1(1,1))
! Total energy
ERHF = ET + EV + EJ + EK
@ -200,7 +234,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
endif
endif
! Level shift
if(level_shift > 0d0 .and. Conv > thresh) then
@ -208,10 +242,20 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
endif
! Diagonalize Fock matrix
Fp = matmul(transpose(X), matmul(F, X))
call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, &
F(1,1), nBas, X(1,1), nBas, &
0.d0, FX(1,1), nBas)
call dgemm("T", "N", nOrb, nOrb, nBas, 1.d0, &
X(1,1), nBas, FX(1,1), nBas, &
0.d0, Fp(1,1), nOrb)
!Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
!c = matmul(X, cp)
call dgemm("N", "N", nBas, nOrb, nOrb, 1.d0, &
X(1,1), nBas, cp(1,1), nOrb, &
0.d0, c(1,1), nBas)
! Density matrix
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
@ -239,6 +283,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
write(*,*)
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
deallocate(tmp1, FX, ERI_chem)
stop
@ -261,5 +306,6 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
end if
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
deallocate(tmp1, FX, ERI_chem)
end subroutine

4
src/HF/RHF_search.f90
View File

@ -126,8 +126,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
ispin = 1
call phLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phRLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phRLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
AB(:,:) = Aph(:,:) + Bph(:,:)

8
src/HF/RHF_stability.f90
View File

@ -38,8 +38,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
ispin = 1
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
call phRLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phRLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
AB(:,:) = A(:,:) + B(:,:)
@ -115,8 +115,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
ispin = 2
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
call phRLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phRLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
AB(:,:) = A(:,:) + B(:,:)

110
src/HF/fermi_dirac_occ.f90 Normal file
View File

@ -0,0 +1,110 @@
subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
! Use Fermi Dirac distribution to set up fractional Occs numbers and adjust the chemical potential to integrate to the N for 2N electron systems
implicit none
! Input variables
integer,intent(in) :: nO,nOrb
double precision,intent(in) :: thrs_N
double precision,intent(in) :: temperature
double precision,intent(in) :: eHF(nOrb)
! Local variables
logical :: backward
integer :: iorb
integer :: isteps
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_1rdm
double precision :: trace_old
! Output variables
double precision,intent(inout):: chem_pot
double precision,intent(inout):: Occ(nOrb)
! Initialize variables
backward = .false.
isteps = 0
delta_chem_pot = 1.0d-1
chem_pot_change = 0d0
grad_electrons = 1d0
trace_1rdm = -1d0
write(*,*)
write(*,*)' Fermi-Dirac distribution for the occ numbers'
write(*,*)
write(*,*)'-------------------------------------'
write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') &
'|','Tr[1D]','|','Chem. Pot.','|'
write(*,*)'-------------------------------------'
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
! First approach close the value with an error lower than 1
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_old=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_old,'|',chem_pot,'|'
do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + delta_chem_pot
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
backward=.true.
chem_pot = chem_pot - delta_chem_pot
delta_chem_pot=-delta_chem_pot
endif
enddo
! Do final search
write(*,*)'-------------------------------------'
isteps=0
delta_chem_pot = 1.0d-1
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + chem_pot_change
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm = sum(Occ(:))
grad_electrons = ( sum(fermi_dirac(eHF,chem_pot+delta_chem_pot,temperature)) &
- sum(fermi_dirac(eHF,chem_pot-delta_chem_pot,temperature)) )/(2.0d0*delta_chem_pot)
chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
enddo
write(*,*)'-------------------------------------'
write(*,*)
write(*,*)
write(*,*) ' Initial occ. numbers '
write(*,*)
do iorb=1,nOrb
write(*,'(3X,F16.10)') Occ(iorb)
enddo
contains
function fermi_dirac(eHF,chem_pot,temperature)
implicit none
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: chem_pot
double precision,intent(in) :: temperature
double precision :: fermi_dirac(nOrb)
fermi_dirac(:) = 1d0 / ( 1d0 + exp((eHF(:) - chem_pot ) / temperature ) )
end function fermi_dirac
end subroutine

159
src/HF/fix_chem_pot.f90 Normal file
View File

@ -0,0 +1,159 @@
subroutine fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_)
! Fix the chemical potential to integrate to the N for 2N electron systems
implicit none
! Input variables
integer,intent(in) :: nO,nOrb,nOrb2,nSCF
double precision,intent(in) :: thrs_N
! Local variables
logical :: backward
integer :: iorb
integer :: isteps
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_2up
double precision :: trace_up
double precision :: trace_down
double precision :: trace_2down
double precision :: trace_old
double precision,allocatable :: R_tmp(:,:)
double precision,allocatable :: cp_tmp(:,:)
! Output variables
double precision :: trace_1rdm
double precision,intent(inout):: chem_pot
double precision,intent(inout):: cp(nOrb2,nOrb2)
double precision,intent(inout):: R(nOrb2,nOrb2)
double precision,intent(inout):: eHFB_(nOrb2)
double precision,intent(inout):: H_hfb(nOrb2,nOrb2)
! Initialize
backward = .false.
isteps = 0
delta_chem_pot = 1.0d-1
chem_pot_change = 0d0
grad_electrons = 1d0
trace_1rdm = -1d0
allocate(R_tmp(nOrb2,nOrb2))
allocate(cp_tmp(nOrb2,nOrb2))
! Set H_HFB to its non-chemical potential dependent contribution
do iorb=1,nOrb
H_hfb(iorb,iorb)=H_hfb(iorb,iorb)+chem_pot
H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)-chem_pot
enddo
write(*,*)
write(*,'(a,I5)') ' Fixing the Tr[1D] at iteration ',nSCF
write(*,*)
write(*,*)'------------------------------------------------------'
write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') &
'|','Tr[1D]','|','Chem. Pot.','|','Grad N','|'
write(*,*)'------------------------------------------------------'
! First approach close the value with an error lower than 1
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_old,H_hfb,cp,R,eHFB_)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
'|',trace_old,'|',chem_pot,'|',grad_electrons,'|'
do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + delta_chem_pot
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|'
if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
backward=.true.
chem_pot = chem_pot - delta_chem_pot
delta_chem_pot=-delta_chem_pot
endif
enddo
! Do final search
write(*,*)'------------------------------------------------------'
isteps = 0
delta_chem_pot = 1.0d-3
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + chem_pot_change
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
call diag_H_hfb(nOrb,nOrb2,chem_pot+2d0*delta_chem_pot,trace_2up,H_hfb,cp_tmp,R_tmp,eHFB_)
call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp_tmp,R_tmp,eHFB_)
call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp_tmp,R_tmp,eHFB_)
call diag_H_hfb(nOrb,nOrb2,chem_pot-2d0*delta_chem_pot,trace_2down,H_hfb,cp_tmp,R_tmp,eHFB_)
! grad_electrons = (trace_up-trace_down)/(2d0*delta_chem_pot)
grad_electrons = (-trace_2up+8d0*trace_up-8d0*trace_down+trace_2down)/(12d0*delta_chem_pot)
chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|'
enddo
write(*,*)'------------------------------------------------------'
write(*,*)
! Reset H_HFB to its chemical potential version
do iorb=1,nOrb
H_hfb(iorb,iorb)=H_hfb(iorb,iorb)-chem_pot
H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)+chem_pot
enddo
deallocate(R_tmp,cp_tmp)
end subroutine
subroutine diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
! Fix the chemical potential to integrate to the N for 2N electron systems
implicit none
! Input variables
integer,intent(in) :: nOrb,nOrb2
double precision,intent(in) :: H_hfb(nOrb2,nOrb2)
double precision,intent(in) :: chem_pot
! Local variables
integer :: iorb
! Output variables
double precision,intent(out) :: trace_1rdm
double precision,intent(inout):: cp(nOrb2,nOrb2)
double precision,intent(inout):: R(nOrb2,nOrb2)
double precision,intent(inout):: eHFB_(nOrb2)
cp(:,:) = H_hfb(:,:)
do iorb=1,nOrb
cp(iorb,iorb) = cp(iorb,iorb) - chem_pot
cp(iorb+nOrb,iorb+nOrb) = cp(iorb+nOrb,iorb+nOrb) + chem_pot
enddo
! Diagonalize H_HFB matrix
call diagonalize_matrix(nOrb2,cp,eHFB_)
! Build R and extract P and Panom
trace_1rdm = 0d0
R(:,:) = 0d0
do iorb=1,nOrb
R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb)))
enddo
do iorb=1,nOrb
trace_1rdm = trace_1rdm + R(iorb,iorb)
enddo
end subroutine

79
src/HF/print_HFB.f90 Normal file
View File

@ -0,0 +1,79 @@
! ---
subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole)
! Print one-electron energies and other stuff
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
double precision,intent(in) :: EJ
double precision,intent(in) :: EK
double precision,intent(in) :: EL
double precision,intent(in) :: ERHF
double precision,intent(in) :: chem_pot
double precision,intent(in) :: N_anom
double precision,intent(in) :: dipole(ncart)
! Local variables
integer :: iorb
integer :: ixyz
double precision :: trace_occ
logical :: dump_orb = .false.
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A33)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au'
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au'
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au'
write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au'
write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au'
write(*,'(A33,1X,F16.10,A3)') ' Anomalous energy = ',EL,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',ERHF,' au'
write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
write(*,'(A33,1X,F16.10,A3)') ' HFB energy = ',ERHF + ENuc,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Chemical potential = ',chem_pot,' au'
write(*,'(A33,1X,F16.10,A3)') ' | Anomalous dens | = ',N_anom,' '
write(*,'(A50)') '---------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '---------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' HFB occupation numbers '
write(*,'(A50)') '---------------------------------------'
trace_occ=0d0
do iorb=1,nOrb
if(abs(Occ(nOrb+1-iorb))>1d-8) then
write(*,'(I7,10F15.8)') iorb,2d0*Occ(nOrb+1-iorb)
endif
trace_occ=trace_occ+2d0*Occ(iorb)
enddo
write(*,*)
write(*,'(A33,1X,F16.10,A3)') ' Trace [ 1D ] = ',trace_occ,' '
write(*,*)
end subroutine

124
src/HF/read_restart_HFB.f90 Normal file
View File

@ -0,0 +1,124 @@
! ---
subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: S(nBas,nBas)
! Local variables
integer :: ibas,iorb,iorb1,iunit=667
integer :: nBas_
integer :: nOrb_
double precision :: chem_pot_
double precision :: max_diff
double precision :: val_read
double precision,allocatable :: eigVAL(:)
double precision,allocatable :: c_tmp(:,:)
double precision,allocatable :: S_mol(:,:)
double precision,allocatable :: X_mol(:,:)
! Output variables
double precision,intent(out) :: chem_pot
double precision,intent(out) :: Occ(nOrb)
double precision,intent(out) :: c(nBas,nOrb)
! Dump results
allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb))
c_tmp=0d0
S_mol=0d0
X_mol=0d0
open(unit=iunit,form='unformatted',file='hfb_bin',status='old')
read(iunit) nBas_,nOrb_
read(iunit) chem_pot_
do iorb=1,nOrb
do ibas=1,nBas
read(iunit) val_read
c_tmp(ibas,iorb) = val_read
enddo
enddo
do iorb=1,nOrb
read(iunit) eigVAL(iorb)
enddo
close(iunit)
if(nBas==nBas_ .and. nOrb==nOrb_) then
write(*,*)
write(*,*)' Reading restart file'
write(*,*)
chem_pot = chem_pot_
Occ(:) = eigVAL(:)
c(:,:) = c_tmp(:,:)
! Check the orthonormality
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff>1d-8) then
write(*,*) ' '
write(*,'(a)') ' Orthonormality violations. Applying Lowdin orthonormalization.'
write(*,*) ' '
X_mol = S_mol
S_mol = 0d0
call diagonalize_matrix(nOrb,X_mol,eigVAL)
do iorb=1,nOrb
S_mol(iorb,iorb) = 1d0/(sqrt(eigVAL(iorb)) + 1d-10)
enddo
X_mol = matmul(X_mol,matmul(S_mol,transpose(X_mol)))
c = matmul(c,X_mol)
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff<=1d-8) then
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
else
write(*,*) ' '
write(*,'(a)') ' Error in Lowdin orthonormalization.'
write(*,*) ' '
stop
endif
else
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
endif
endif
deallocate(eigVAL,c_tmp,S_mol,X_mol)
end subroutine

90
src/HF/reshufle_hfb.f90 Normal file
View File

@ -0,0 +1,90 @@
! ---
subroutine reshufle_hfb(nBas2,nOrb2,c_hfb,eHFB,project)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas2, nOrb2
double precision,intent(inout) :: eHFB(nOrb2)
double precision,intent(inout) :: project(nOrb2)
double precision,intent(inout) :: c_hfb(nBas2,nOrb2)
! Local variables
integer :: iorb,iorb2,ipos
double precision :: val
double precision,allocatable :: eHFB_tmp(:)
double precision,allocatable :: project_tmp(:)
double precision,allocatable :: c_tmp(:,:)
allocate(c_tmp(nBas2,nOrb2),project_tmp(nOrb2),eHFB_tmp(nOrb2))
!write(*,*) 'e_HFB'
!do iorb=1,nOrb2
!write(*,*) eHFB(iorb),project(iorb)
!enddo
! Reshufle using MOM
do iorb=1,nOrb2
val=-1d0
ipos=0
do iorb2=1,nOrb2
if(project(iorb2)>val) then
val=project(iorb2)
ipos=iorb2
endif
enddo
project(ipos)=-1.0d1
project_tmp(iorb)=val
eHFB_tmp(iorb)=eHFB(ipos)
c_tmp(:,iorb)=c_hfb(:,ipos)
enddo
! Reshufle using energies
do iorb=1,nOrb2/2 ! electronic 1-RDM
val=1d10
ipos=0
do iorb2=1,nOrb2/2
if(eHFB_tmp(iorb2)<val) then
val=eHFB_tmp(iorb2)
ipos=iorb2
endif
enddo
eHFB_tmp(ipos)=2d10
eHFB(iorb)=val
project(iorb)=project_tmp(ipos)
c_hfb(:,iorb)=c_tmp(:,ipos)
enddo
do iorb=nOrb2/2+1,nOrb2 ! hole 1-RDM
val=1d10
ipos=0
do iorb2=nOrb2/2+1,nOrb2
if(eHFB_tmp(iorb2)<val) then
val=eHFB_tmp(iorb2)
ipos=iorb2
endif
enddo
eHFB_tmp(ipos)=2d10
eHFB(iorb)=val
project(iorb)=project_tmp(ipos)
c_hfb(:,iorb)=c_tmp(:,ipos)
enddo
!write(*,*) 'e_HFB'
!do iorb=1,nOrb2/2
!write(*,*) eHFB(iorb),project(iorb)
!enddo
!write(*,*) '------'
!do iorb=nOrb2/2+1,nOrb2
!write(*,*) eHFB(iorb),project(iorb)
!enddo
deallocate(c_tmp,project_tmp,eHFB_tmp)
end subroutine

44
src/HF/write_restart_HFB.f90 Normal file
View File

@ -0,0 +1,44 @@
! ---
subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: chem_pot
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
! Local variables
integer :: ibas,iorb,iunit=666
double precision :: val_print
! Dump results
open(unit=iunit,form='unformatted',file='hfb_bin')
write(iunit) nBas,nOrb
write(iunit) chem_pot
do iorb=1,nOrb
do ibas=1,nBas
val_print = c(ibas,nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo
enddo
do iorb=1,nOrb
val_print = Occ(nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo
write(iunit) iunit
close(iunit)
end subroutine

2
src/LR/phLR.f90 → src/LR/phRLR.f90
View File

@ -1,4 +1,4 @@
subroutine phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute linear response

2
src/LR/phLR_A.f90 → src/LR/phRLR_A.f90
View File

@ -1,4 +1,4 @@
subroutine phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
subroutine phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
! Compute resonant block of the ph channel

2
src/LR/phLR_B.f90 → src/LR/phRLR_B.f90
View File

@ -1,4 +1,4 @@
subroutine phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
subroutine phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
! Compute the coupling block of the ph channel

61
src/LR/ppGLR.f90
View File

@ -1,7 +1,7 @@
subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
!
! Solve the pp-RPA linear eigenvalue problem
! Solve the pp-RPA linear eigenvalue problem for a generalized reference
!
! right eigen-problem: H R = R w
! left eigen-problem: H.T L = L w
@ -49,6 +49,7 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP))
! Hermitian case for TDA
if(TDA) then
X1(:,:) = +Cpp(:,:)
@ -71,7 +72,7 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
M( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO)
M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO))
! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! Diagonalize the pp matrix
@ -81,36 +82,36 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! else
! else
! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
! imp_bio = .True. ! impose bi-orthogonality
! verbose = .False.
! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
! imp_bio = .True. ! impose bi-orthogonality
! verbose = .False.
! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! endif
! endif
end if

4
src/LR/ppLR_transition_vectors.f90
View File

@ -27,8 +27,8 @@ subroutine ppLR_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_
integer :: a,b,c,d,ab,cd
integer :: i,j,k,l,ij,kl
integer :: maxOO = 10
integer :: maxVV = 10
integer :: maxOO = 1000
integer :: maxVV = 0
double precision :: S2
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: os1(:)

2
src/LR/ppLR.f90 → src/LR/ppRLR.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
subroutine ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
!
! Solve the pp-RPA linear eigenvalue problem

2
src/LR/ppLR_B.f90 → src/LR/ppRLR_B.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
subroutine ppRLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
! Compute the B matrix of the pp channel

2
src/LR/ppLR_C.f90 → src/LR/ppRLR_C.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
subroutine ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
! Compute the C matrix of the pp channel

2
src/LR/ppLR_C_od.f90 → src/LR/ppRLR_C_od.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Cpp)
subroutine ppRLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Cpp)
! Compute the C matrix of the pp channel (without diagonal term)

2
src/LR/ppLR_D.f90 → src/LR/ppRLR_D.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
subroutine ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
! Compute the D matrix of the pp channel

2
src/LR/ppLR_D_od.f90 → src/LR/ppRLR_D_od.f90
View File

@ -1,4 +1,4 @@
subroutine ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Dpp)
subroutine ppRLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Dpp)
! Compute the D matrix of the pp channel (without the diagonal term)

118
src/QuAcK/BQuAcK.f90 Normal file
View File

@ -0,0 +1,118 @@
subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,temperature,sigma,chem_pot_hf,restart_hfb)
! Restricted branch of QuAcK
implicit none
include 'parameters.h'
character(len=256),intent(in) :: working_dir
logical,intent(in) :: dotest
logical,intent(in) :: doHFB
logical,intent(in) :: restart_hfb
logical,intent(in) :: chem_pot_hf
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature,sigma
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
integer,intent(in) :: guess_type
! Local variables
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:)
double precision,allocatable :: PanomHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision,allocatable :: Delta(:,:)
double precision :: ERHF,EHFB
double precision,allocatable :: ERI_AO(:,:,:,:)
write(*,*)
write(*,*) '******************************'
write(*,*) '* Bogoliubov Branch of QuAcK *'
write(*,*) '******************************'
write(*,*)
!-------------------!
! Memory allocation !
!-------------------!
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas))
allocate(PanomHF(nBas,nBas))
allocate(FHF(nBas,nBas))
allocate(Delta(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!--------------------------------!
! Hartree-Fock Bogoliubov module !
!--------------------------------!
if(doHFB) then
! Run first a RHF calculation
call wall_time(start_HF)
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
! Continue with a HFB calculation
call wall_time(start_HF)
call HFB(dotest,maxSCF_HF,thresh_HF,max_diis_HF,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,PanomHF, &
FHF,Delta,temperature,sigma,chem_pot_hf,restart_hfb)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HFB = ',t_HF,' seconds'
write(*,*)
end if
! Memory deallocation
deallocate(eHF)
deallocate(cHF)
deallocate(PHF)
deallocate(PanomHF)
deallocate(FHF)
deallocate(Delta)
deallocate(ERI_AO)
end subroutine

34
src/QuAcK/QuAcK.f90
View File

@ -3,8 +3,8 @@ program QuAcK
implicit none
include 'parameters.h'
logical :: doRQuAcK,doUQuAcK,doGQuAcK
logical :: doRHF,doUHF,doGHF,doROHF
logical :: doRQuAcK,doUQuAcK,doGQuAcK,doBQuAcK
logical :: doRHF,doUHF,doGHF,doROHF,doHFB
logical :: dostab,dosearch
logical :: doMP2,doMP3
logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
@ -74,6 +74,10 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
logical :: chem_pot_hf
logical :: restart_hfb
double precision :: temperature,sigma
character(len=256) :: working_dir
! Check if the right number of arguments is provided
@ -109,7 +113,7 @@ program QuAcK
!------------------!
call read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doRHF,doUHF,doGHF,doROHF,doHFB, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD, &
@ -136,7 +140,8 @@ program QuAcK
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf,restart_hfb)
!------------------!
! Hardware !
@ -214,6 +219,9 @@ program QuAcK
doGQuAcK = .false.
if(doGHF) doGQuAcK = .true.
doBQuAcK = .false.
if(doHFB) doBQuAcK = .true.
!-----------------!
! Initialize Test !
!-----------------!
@ -283,6 +291,15 @@ program QuAcK
maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
!--------------------------!
! Bogoliubov QuAcK branch !
!--------------------------!
if(doBQuAcK) &
call BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix, &
temperature,sigma,chem_pot_hf,restart_hfb)
!-----------!
! Stop Test !
!-----------!
@ -305,10 +322,9 @@ program QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds'
write(*,*)
deallocate(S)
deallocate(T)
deallocate(V)
deallocate(Hc)
deallocate(dipole_int_AO)
! Memory deallocation
deallocate(ZNuc,rNuc)
deallocate(S,T,V,Hc,dipole_int_AO)
end program

1
src/QuAcK/RQuAcK.f90
View File

@ -380,6 +380,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,
end if
! Memory deallocation
deallocate(eHF)
deallocate(cHF)
deallocate(PHF)

8
src/QuAcK/read_methods.f90
View File

@ -1,5 +1,5 @@
subroutine read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doRHF,doUHF,doGHF,doROHF,doHFB, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
@ -23,7 +23,7 @@ subroutine read_methods(working_dir, &
! Output variables
logical,intent(out) :: doRHF,doUHF,doGHF,doROHF
logical,intent(out) :: doRHF,doUHF,doGHF,doROHF,doHFB
logical,intent(out) :: doMP2,doMP3
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
@ -60,13 +60,15 @@ subroutine read_methods(working_dir, &
doUHF = .false.
doGHF = .false.
doROHF = .false.
doHFB = .false.
read(1,*)
read(1,*) ans1,ans2,ans3,ans4
read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') doRHF = .true.
if(ans2 == 'T') doUHF = .true.
if(ans3 == 'T') doGHF = .true.
if(ans4 == 'T') doROHF = .true.
if(ans5 == 'T') doHFB = .true.
! Read MPn methods

21
src/QuAcK/read_options.f90
View File

@ -7,7 +7,8 @@ subroutine read_options(working_dir,
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf,restart_hfb)
! Read desired methods
@ -72,6 +73,11 @@ subroutine read_options(working_dir,
logical,intent(out) :: dBSE
logical,intent(out) :: dTDA
logical,intent(out) :: chem_pot_hf
logical,intent(out) :: restart_hfb
double precision,intent(out) :: temperature
double precision,intent(out) :: sigma
! Local variables
character(len=1) :: ans1,ans2,ans3,ans4,ans5
@ -217,6 +223,19 @@ subroutine read_options(working_dir,
if(ans4 == 'T') dBSE = .true.
if(ans5 == 'F') dTDA = .false.
! Options for Hartree-Fock Bogoliubov
temperature = 0d0
sigma = 1d0
chem_pot_hf = .false.
restart_hfb = .false.
read(1,*)
read(1,*) temperature,sigma,ans1,ans2
if(ans1 == 'T') chem_pot_hf = .true.
if(ans2 == 'T') restart_hfb = .true.
endif
! Close file with options

12
src/RPA/crACFDT.f90
View File

@ -87,10 +87,10 @@ subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
lambda = -rAC(iAC)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
@ -128,10 +128,10 @@ subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
lambda = -rAC(iAC)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))

4
src/RPA/crGRPA.f90
View File

@ -58,8 +58,8 @@ subroutine crGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
call phLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call print_excitation_energies('crRPA@GHF','spinorbital',nS,Om)

12
src/RPA/crRRPA.f90
View File

@ -69,10 +69,10 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('crRPA@RHF','singlet',nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -84,10 +84,10 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('crRPA@RHF','triplet',nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)

12
src/RPA/phACFDT.f90
View File

@ -97,10 +97,10 @@ subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
lambda = rAC(iAC)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
@ -138,10 +138,10 @@ subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
lambda = rAC(iAC)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))

2
src/RPA/phGRPAx.f90
View File

@ -59,7 +59,7 @@ subroutine phGRPAx(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call print_excitation_energies('phRPAx@GHF','spinorbital',nS,Om)

19
src/RPA/phRRPA.f90
View File

@ -73,17 +73,10 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 1
!call wall_time(t1)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
!call wall_time(t2)
!print *, "wall time for dRPA on CPU (sec) = ", t2 - t1
!do ia = 1, nS
! write(112, *) Om(ia)
!enddo
!stop
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -95,10 +88,10 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)

14
src/RPA/phRRPAx.f90
View File

@ -70,10 +70,10 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPAx@RHF','singlet',nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -85,10 +85,10 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPAx@RHF','triplet',nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -99,7 +99,7 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
end if
endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

16
src/RPA/ppACFDT.f90
View File

@ -89,11 +89,11 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC)
lambda = rAC(iAC)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))
@ -138,11 +138,11 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC)
lambda = rAC(iAC)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))

6
src/RPA/ppACFDT_correlation_energy.f90
View File

@ -39,9 +39,9 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X
! Build pp matrices
call ppLR_B (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B)
call ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C)
call ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D)
call ppRLR_B (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B)
call ppRLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C)
call ppRLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D)
! Compute Tr(K x P_lambda)

16
src/RPA/ppRRPA.f90
View File

@ -70,11 +70,11 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
@ -102,11 +102,11 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
if(.not.TDA) call ppRLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppRLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppRLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
call ppRLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)

32
src/utils/sort.f90
View File

@ -60,6 +60,38 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level)
end if
end subroutine
subroutine set_order_LR(x,y,iorder,isize,jsize)
implicit none
integer :: isize,jsize
double precision :: x(isize,jsize)
double precision :: y(isize,jsize)
double precision,allocatable :: xtmp(:,:)
double precision,allocatable :: ytmp(:,:)
integer :: iorder(*)
integer :: i,j
allocate(xtmp(isize,jsize),ytmp(isize,jsize))
do i=1,isize
do j=1,jsize
xtmp(i,j) = x(i,iorder(j))
ytmp(i,j) = y(i,iorder(j))
end do
end do
do i=1,isize
do j=1,jsize
x(i,j) = xtmp(i,j)
y(i,j) = ytmp(i,j)
end do
end do
deallocate(xtmp,ytmp)
end subroutine
subroutine set_order(x,iorder,isize,jsize)
implicit none

6
src/utils/utils.f90
View File

@ -929,9 +929,11 @@ integer*8 function Yoshimine_2ind(a, b)
integer*8, intent(in) :: a, b
if(a > b) then
Yoshimine_2ind = (a * (a - 1)) / 2 + b
!Yoshimine_2ind = (a * (a - 1)) / 2 + b
Yoshimine_2ind = shiftr(a * (a - 1), 1) + b
else
Yoshimine_2ind = (b * (b - 1)) / 2 + a
!Yoshimine_2ind = (b * (b - 1)) / 2 + a
Yoshimine_2ind = shiftr(b * (b - 1), 1) + a
endif
return

52
src/utils/wrap_lapack.f90
View File

@ -1,3 +1,53 @@
subroutine diagonalize_general_matrix_LR(N,A,WR,VL,VR)
! Diagonalize a non-symmetric square matrix
implicit none
! Input variables
integer,intent(in) :: N
double precision,intent(inout):: A(N,N)
double precision,intent(out) :: VL(N,N)
double precision,intent(out) :: VR(N,N)
double precision,intent(out) :: WR(N)
! Local variables
integer :: i
double precision :: tmp
integer :: lwork,info
double precision,allocatable :: work(:),WI(:)
! Memory allocation
allocate(work(1),WI(N))
lwork = -1
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
do i=1,N
tmp = dot_product(vl(:,i),vr(:,i))
vl(:,i) = vl(:,i)/tmp
end do
call matout(N,N,matmul(transpose(VL),VR))
deallocate(work,WI)
if(info /= 0) then
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
end if
end subroutine
subroutine diagonalize_general_matrix(N,A,WR,VR)
! Diagonalize a non-symmetric square matrix
@ -31,7 +81,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
deallocate(work, WI, VL)
deallocate(work,WI,VL)
if(info /= 0) then
print*,'Problem in diagonalize_general_matrix (dgeev)!!'