PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-03 18:06:03 +02:00
Code Releases Activity

debug pCCD

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-22 15:20:42 +01:00
parent eee1a1c6b4
commit f43cf3f4e4
7 changed files with 131 additions and 112 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
input/basis
View File

@ -1,9 +1,30 @@
1 3 1 6
S 3 S 8
1 38.3600000 0.0238090 1 17880.0000000 0.0007380
2 5.7700000 0.1548910 2 2683.0000000 0.0056770
3 1.2400000 0.4699870 3 611.5000000 0.0288830
4 173.5000000 0.1085400
5 56.6400000 0.2909070
6 20.4200000 0.4483240
7 7.8100000 0.2580260
8 1.6530000 0.0150630
S 8
1 17880.0000000 -0.0001720
2 2683.0000000 -0.0013570
3 611.5000000 -0.0067370
4 173.5000000 -0.0276630
5 56.6400000 -0.0762080
6 20.4200000 -0.1752270
7 7.8100000 -0.1070380
8 1.6530000 0.5670500
S 1 S 1
1 0.2976000 1.0000000 1 0.4869000 1.0000000
P 3
1 28.3900000 0.0460870
2 6.2700000 0.2401810
3 1.6950000 0.5087440
P 1 P 1
1 1.2750000 1.0000000 1 0.4317000 1.0000000
D 1
1 2.2020000 1.0000000

4
input/methods
View File

@ -3,9 +3,9 @@
# MP2 MP3 MP2-F12 # MP2 MP3 MP2-F12
T F F T F F
# CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD # CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD
F F F T T T T F F F F F F T
# CIS RPA RPAx ppRPA ADC # CIS RPA RPAx ppRPA ADC
F T T T F F F F F F
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3
F F F F F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW

4
input/molecule
View File

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

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1 1
He 0.0000000000 0.0000000000 0.0000000000 Ne 0.0000000000 0.0000000000 0.0000000000

35
input/weight
View File

@ -1,9 +1,30 @@
1 3 1 6
S 3 S 8
1 38.3600000 0.0238090 1 17880.0000000 0.0007380
2 5.7700000 0.1548910 2 2683.0000000 0.0056770
3 1.2400000 0.4699870 3 611.5000000 0.0288830
4 173.5000000 0.1085400
5 56.6400000 0.2909070
6 20.4200000 0.4483240
7 7.8100000 0.2580260
8 1.6530000 0.0150630
S 8
1 17880.0000000 -0.0001720
2 2683.0000000 -0.0013570
3 611.5000000 -0.0067370
4 173.5000000 -0.0276630
5 56.6400000 -0.0762080
6 20.4200000 -0.1752270
7 7.8100000 -0.1070380
8 1.6530000 0.5670500
S 1 S 1
1 0.2976000 1.0000000 1 0.4869000 1.0000000
P 3
1 28.3900000 0.0460870
2 6.2700000 0.2401810
3 1.6950000 0.5087440
P 1 P 1
1 1.2750000 1.0000000 1 0.4317000 1.0000000
D 1
1 2.2020000 1.0000000

2
src/QuAcK/QuAcK.f90
View File

@ -447,7 +447,7 @@ program QuAcK
if(do_pCCD) then if(do_pCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nO(1),nV(1),ERI_MO_basis,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD t_CCD = end_CCD - start_CCD

161
src/QuAcK/pCCD.f90
View File

@ -1,4 +1,4 @@
subroutine pCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
! pair CCD module ! pair CCD module
@ -10,110 +10,86 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas,nEl integer,intent(in) :: nBas,nO,nV
double precision,intent(in) :: ENuc,ERHF double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
integer :: nBas2 integer :: i,j,a,b
integer :: nO
integer :: nV
integer :: nSCF integer :: nSCF
double precision :: Conv double precision :: Conv
double precision :: EcMP2,EcMP3,EcMP4
double precision :: ECCD,EcCCD double precision :: ECCD,EcCCD
double precision,allocatable :: seHF(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: eO(:) double precision,allocatable :: delta_OOVV(:,:)
double precision,allocatable :: eV(:)
double precision,allocatable :: delta_OOVV(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:) double precision,allocatable :: OOOO(:,:)
double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OOVV(:,:)
double precision,allocatable :: OVOV(:,:,:,:) double precision,allocatable :: OVOV(:,:)
double precision,allocatable :: VVVV(:,:,:,:) double precision,allocatable :: OVVO(:,:)
double precision,allocatable :: VVVV(:,:)
double precision,allocatable :: X1(:,:,:,:) double precision,allocatable :: X(:,:)
double precision,allocatable :: X2(:,:) double precision,allocatable :: Y(:,:)
double precision,allocatable :: X3(:,:)
double precision,allocatable :: X4(:,:,:,:)
double precision,allocatable :: u(:,:,:,:) double precision,allocatable :: r(:,:)
double precision,allocatable :: v(:,:,:,:) double precision,allocatable :: t(:,:)
double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:)
double precision,external :: trace_matrix
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'**************************************' write(*,*)'**************************************'
write(*,*)'| CCD calculation |' write(*,*)'| pair CCD calculation |'
write(*,*)'**************************************' write(*,*)'**************************************'
write(*,*) write(*,*)
! Spatial to spin orbitals
nBas2 = 2*nBas
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
deallocate(sERI)
! Define occupied and virtual spaces
nO = 2*nEl
nV = nBas2 - nO
! Form energy denominator ! Form energy denominator
allocate(eO(nO),eV(nV)) allocate(delta_OOVV(nO,nV))
allocate(delta_OOVV(nO,nO,nV,nV))
eO(:) = seHF(1:nO) do i=1,nO
eV(:) = seHF(nO+1:nBas2) do a=1,nV
delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i))
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) enddo
enddo
deallocate(seHF)
! Create integral batches ! Create integral batches
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV)) allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV))
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) do i=1,nO
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2) do j=1,nO
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2) OOOO(i,j) = ERI(i,i,j,j)
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) end do
end do
do i=1,nO
do a=1,nV
OOVV(i,a) = ERI(i,i,nO+a,nO+a)
OVOV(i,a) = ERI(i,nO+a,i,nO+a)
OVVO(i,a) = ERI(i,nO+a,nO+a,i)
end do
end do
do a=1,nV
do b=1,nV
VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b)
end do
end do
deallocate(dbERI)
! MP2 guess amplitudes ! MP2 guess amplitudes
allocate(t2(nO,nO,nV,nV)) allocate(t(nO,nV))
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) t(:,:) = - OOVV(:,:)/delta_OOVV(:,:)
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
EcMP4 = 0d0
! Initialization ! Initialization
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) allocate(r(nO,nV),X(nV,nV),Y(nO,nO))
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
Conv = 1d0 Conv = 1d0
nSCF = 0 nSCF = 0
@ -131,39 +107,40 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
do while(Conv > thresh .and. nSCF < maxSCF) do while(Conv > thresh .and. nSCF < maxSCF)
! Increment ! Increment
nSCF = nSCF + 1 nSCF = nSCF + 1
! Form linear array ! Form intermediate array
X(:,:) = matmul(transpose(OOVV(:,:)),t(:,:))
Y(:,:) = matmul(t(:,:),transpose(OOVV(:,:)))
! Compute residual
call form_u(nO,nV,OOOO,VVVV,OVOV,t2,u) do i=1,nO
do a=1,nV
r(i,a) = - 2d0*(X(a,a) + Y(i,i))*t(i,a)
end do
end do
! Form interemediate arrays r(:,:) = r(:,:) + OOVV(:,:) + delta_OOVV(:,:)*t(:,:) &
- 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t(:,:))*t(:,:) &
+ matmul(t(:,:),transpose(VVVV(:,:))) &
+ matmul(transpose(OOOO(:,:)),t(:,:)) &
+ matmul(Y(:,:),t)
call form_X(nO,nV,OOVV,t2,X1,X2,X3,X4) ! Check convergence
! Form quadratic array Conv = maxval(abs(r(:,:)))
call form_v(nO,nV,X1,X2,X3,X4,t2,v)
! Compute residual
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:)
! Check convergence
Conv = maxval(abs(r2(:,:,:,:)))
! Update amplitudes ! Update amplitudes
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) t(:,:) = t(:,:) - r(:,:)/delta_OOVV(:,:)
! Compute correlation energy ! Compute correlation energy
EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) EcCCD = trace_matrix(nO,matmul(t(:,:),transpose(OOVV(:,:))))
if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.))
! Dump results ! Dump results