PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-23 19:27:38 +02:00
Code Releases Activity

BSE clean up

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-11 16:41:20 +01:00
parent 08d97376f2
commit 63e6be3b2c
22 changed files with 214 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.CO
View File

@ -2,4 +2,4 @@
2 7 7 0 0
# Znuc x y z
C 0. 0. 0.
O 0. 0. 2.106
O 0. 0. 2.134

2
examples/molecule.F2
View File

@ -2,4 +2,4 @@
2 9 9 0 0
# Znuc x y z
F 0. 0. 0.
F 0. 0. 2.6682933
F 0. 0. 2.622798396593

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
H 0. 0. 1.399

2
examples/molecule.LiH
View File

@ -2,4 +2,4 @@
2 2 2 0 0
# Znuc x y z
Li 0. 0. 0.
H 0. 0. 3.0141132
H 0. 0. 3.017

97
input/basis
View File

@ -1,49 +1,72 @@
1 9
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
1 15
S 9
1 6601.0000000 0.0001170
2 989.7000000 0.0009110
3 225.7000000 0.0047280
4 64.2900000 0.0191970
5 21.1800000 0.0630470
6 7.7240000 0.1632080
7 3.0030000 0.3148270
8 1.2120000 0.3939360
9 0.4930000 0.1969180
S 9
1 6601.0000000 -0.0000180
2 989.7000000 -0.0001420
3 225.7000000 -0.0007410
4 64.2900000 -0.0030200
5 21.1800000 -0.0101230
6 7.7240000 -0.0270940
7 3.0030000 -0.0573590
8 1.2120000 -0.0938950
9 0.4930000 -0.1210910
S 1
1 0.0280500 1.0000000
1 0.0951500 1.0000000
S 1
1 0.0086400 1.0000000
1 0.0479100 1.0000000
S 1
1 0.0222000 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
1 6.2500000 0.0033880
2 1.3700000 0.0193160
3 0.3672000 0.0791040
P 1
1 0.0240300 1.0000000
1 0.1192000 1.0000000
P 1
1 0.0057900 1.0000000
1 0.0447400 1.0000000
P 1
1 0.0179500 1.0000000
D 1
1 0.1239000 1.0000000
1 0.3440000 1.0000000
D 1
1 0.0725000 1.0000000
2 5
1 0.1530000 1.0000000
D 1
1 0.0680000 1.0000000
F 1
1 0.2460000 1.0000000
F 1
1 0.1292000 1.0000000
G 1
1 0.2380000 1.0000000
2 10
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
F 1
1 1.3970000 1.0000000

6
input/methods
View File

@ -7,10 +7,10 @@
# CIS RPA RPAx ppRPA ADC
F F F F F
# GF2 GF3
T F
F F
# G0W0 evGW qsGW
F F F
# G0T0 evGT qsGT
T F F
# G0T0 evGT qsGT
F F F
# MCMP2
F

2
input/molecule
View File

@ -2,4 +2,4 @@
2 2 2 0 0
# Znuc x y z
Li 0. 0. 0.
H 0. 0. 3.0141132
H 0. 0. 3.017

2
input/molecule.xyz
View File

@ -1,4 +1,4 @@
2
Li 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 1.5950001314
H 0.0000000000 0.0000000000 1.5965277602

4
input/options
View File

@ -7,10 +7,10 @@
# CIS/TDHF/BSE: singlet triplet
T T
# GF: maxSCF thresh DIIS n_diis renormalization
64 0.00001 T 5 3
1 0.00001 T 5 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F F 0.000
# ACFDT: AC Kx XBS
T F T
T F F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

97
input/weight
View File

@ -1,49 +1,72 @@
1 9
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
1 15
S 9
1 6601.0000000 0.0001170
2 989.7000000 0.0009110
3 225.7000000 0.0047280
4 64.2900000 0.0191970
5 21.1800000 0.0630470
6 7.7240000 0.1632080
7 3.0030000 0.3148270
8 1.2120000 0.3939360
9 0.4930000 0.1969180
S 9
1 6601.0000000 -0.0000180
2 989.7000000 -0.0001420
3 225.7000000 -0.0007410
4 64.2900000 -0.0030200
5 21.1800000 -0.0101230
6 7.7240000 -0.0270940
7 3.0030000 -0.0573590
8 1.2120000 -0.0938950
9 0.4930000 -0.1210910
S 1
1 0.0280500 1.0000000
1 0.0951500 1.0000000
S 1
1 0.0086400 1.0000000
1 0.0479100 1.0000000
S 1
1 0.0222000 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
1 6.2500000 0.0033880
2 1.3700000 0.0193160
3 0.3672000 0.0791040
P 1
1 0.0240300 1.0000000
1 0.1192000 1.0000000
P 1
1 0.0057900 1.0000000
1 0.0447400 1.0000000
P 1
1 0.0179500 1.0000000
D 1
1 0.1239000 1.0000000
1 0.3440000 1.0000000
D 1
1 0.0725000 1.0000000
2 5
1 0.1530000 1.0000000
D 1
1 0.0680000 1.0000000
F 1
1 0.2460000 1.0000000
F 1
1 0.1292000 1.0000000
G 1
1 0.2380000 1.0000000
2 10
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
F 1
1 1.3970000 1.0000000

15
src/QuAcK/ACFDT.f90
View File

@ -1,5 +1,5 @@
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eHF,e,Omega,XpY,XmY,rho,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -19,6 +19,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
@ -76,16 +77,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
if(doXBS) then
call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -121,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
if(doXBS) then
call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
@ -130,7 +132,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)

5
src/QuAcK/Bethe_Salpeter_A_matrix.f90
View File

@ -39,6 +39,7 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
enddo
! A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi
enddo
@ -46,7 +47,7 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
enddo
enddo
print*,'BSE A'
call matout(nS,nS,A_lr)
! print*,'BSE A'
! call matout(nS,nS,A_lr)
end subroutine Bethe_Salpeter_A_matrix

5
src/QuAcK/Bethe_Salpeter_B_matrix.f90
View File

@ -39,6 +39,7 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
enddo
! B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi
enddo
@ -46,7 +47,7 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
enddo
enddo
! print*,'BSE B'
! call matout(nS,nS,B_lr)
print*,'BSE B'
call matout(nS,nS,B_lr)
end subroutine Bethe_Salpeter_B_matrix

11
src/QuAcK/G0T0.f90
View File

@ -1,4 +1,4 @@
subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF,eG0T0)
subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Perform G0W0 calculation with a T-matrix self-energy (G0T0)
@ -7,9 +7,6 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
! Input variables
logical,intent(in) :: BSE
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
@ -36,9 +33,9 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
double precision,allocatable :: SigT(:)
double precision,allocatable :: Z(:)
! Output variables
double precision,allocatable :: eG0T0(:)
double precision,intent(out) :: eG0T0(nBas)
! Output variables
! Hello world
@ -64,7 +61,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nO-nC,nVVt),rho2t(nBas,nV-nR,nOOt), &
SigT(nBas),Z(nBas))
SigT(nBas),Z(nBas),eG0T0(nBas))
!----------------------------------------------
! Singlet manifold

2
src/QuAcK/G0W0.f90
View File

@ -166,7 +166,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC)
if(exchange_kernel) then

3
src/QuAcK/QuAcK.f90
View File

@ -597,8 +597,7 @@ program QuAcK
if(doG0T0) then
call cpu_time(start_G0T0)
! call G0T0(BSE,singlet_manifold,triplet_manifold,eta, &
! nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
! call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_G0T0)

2
src/QuAcK/evGW.f90
View File

@ -257,7 +257,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

6
src/QuAcK/linear_response_A_matrix.f90
View File

@ -47,8 +47,10 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
jb = jb + 1
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
! + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
! - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+ (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*lambda*ERI(i,a,j,b)
enddo
enddo

6
src/QuAcK/linear_response_B_matrix.f90
View File

@ -44,8 +44,10 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
do b=nO+1,nBas-nR
jb = jb + 1
B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
! B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
! - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*lambda*ERI(i,a,b,j)
enddo
enddo

64
src/QuAcK/linear_response_pp.f90
View File

@ -17,6 +17,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
! Local variables
integer :: ab,cd,ij,kl
integer :: p,q,r,s
double precision :: trace_matrix
double precision,allocatable :: B(:,:)
double precision,allocatable :: C(:,:)
@ -64,31 +65,46 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO))
open(unit=42,file='B.dat')
open(unit=43,file='C.dat')
open(unit=44,file='D.dat')
open(unit=42,file='B.dat')
open(unit=43,file='C.dat')
open(unit=44,file='D.dat')
open(unit=45,file='ERI.dat')
open(unit=46,file='eps.dat')
do ab=1,nVV
do ij=1,nOO
if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij)
end do
end do
do ab=1,nVV
do cd=1,nVV
if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd)
end do
end do
do ij=1,nOO
do kl=1,nOO
if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl)
end do
end do
do p=1,nBas
write(46,*) p,e(p)
do q=1,nBas
do r=1,nBas
do s=1,nBas
if(abs(ERI(p,q,r,s)) > 1d-15) write(45,*) p,q,r,s,ERI(p,q,r,s)
end do
end do
end do
end do
do ab=1,nVV
do ij=1,nOO
if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij)
end do
end do
do ab=1,nVV
do cd=1,nVV
if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd)
end do
end do
do ij=1,nOO
do kl=1,nOO
if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl)
end do
end do
close(42)
close(43)
close(44)
close(42)
close(43)
close(44)
close(45)
close(46)
! print*, 'pp-RPA matrix'

2
src/QuAcK/qsGW.f90
View File

@ -288,7 +288,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho,EcAC)
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

19
src/QuAcK/sort_ppRPA.f90
View File

@ -101,7 +101,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Orthogonalize eigenvectors
S1 = matmul(transpose(Z1),matmul(M,Z1))
S1 = + matmul(transpose(Z1),matmul(M,Z1))
S2 = - matmul(transpose(Z2),matmul(M,Z2))
call orthogonalization_matrix(1,nVV,S1,O1)
@ -116,12 +116,17 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
X2(1:nVV,1:nOO) = Z2( 1: nVV,1:nOO)
Y2(1:nOO,1:nOO) = Z2(nVV+1:nOO+nVV,1:nOO)
write(*,*) 'X1t.X1 - Y1t.Y1'
call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
write(*,*) 'X2t.X2 - Y2t.Y2'
call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
write(*,*) 'X1t.X2 - Y1t.Y2'
call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
! write(*,*) 'Z1t.M.Z1'
! call matout(nVV,nVV,matmul(matmul(transpose(Z1),M),Z1))
! write(*,*) 'Z2t.M.Z2'
! call matout(nOO,nOO,matmul(matmul(transpose(Z2),M),Z2))
! write(*,*) 'X1t.X1 - Y1t.Y1'
! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
! write(*,*) 'X2t.X2 - Y2t.Y2'
! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
! write(*,*) 'X1t.X2 - Y1t.Y2'
! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
end subroutine sort_ppRPA