PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-04 05:03:49 +01:00
Code Releases Activity

removing SP files

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-03-19 10:15:55 +01:00
parent 4c821fe4ac
commit a8dba8205c
7 changed files with 36 additions and 509 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
src/QuAcK/Makefile Normal file
View File

@ -0,0 +1,36 @@
IDIR =../../include
BDIR =../../bin
ODIR = obj
OODIR = ../IntPak/obj
SDIR =.
FC = gfortran -I$(IDIR)
ifeq ($(DEBUG),1)
FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
else
FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
SRCF90 = $(wildcard *.f90)
SRC = $(wildcard *.f)
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
$(ODIR)/%.o: %.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.f
$(FC) -c -o $@ $< $(FFLAGS)
$(BDIR)/QuAcK: $(OBJ)
$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
debug:
DEBUG=1 make $(BDIR)/QuAcK
#DEBUG=1 make clean $(BDIR)/QuAcK
clean:
rm -f $(ODIR)/*.o $(BDIR)/QuAcK $(BDIR)/debug

170
src/QuAcK/SPHF.f90
View File

@ -1,170 +0,0 @@
subroutine SPHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Perform restricted Hartree-Fock calculation
implicit none
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
double precision,intent(in) :: thresh
integer,intent(in) :: nBas,nO
double precision,intent(in) :: ENuc
double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas)
! Local variables
integer :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
! Output variables
double precision,intent(out) :: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Restricted Hartree-Fock calculation |'
write(*,*)'************************************************'
write(*,*)
! Useful quantities
nBasSq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
elseif(guess_type == 2) then
call random_number(c)
endif
P(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| SPHF calculation |'
write(*,*)'----------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','HF energy','|','Conv','|','HL Gap','|'
write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock matrix
call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 2d0*K(:,:)
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! Density matrix
P(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap
if(nBas > nO) then
Gap = e(nO+1) - e(nO)
else
Gap = 0d0
endif
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|'
enddo
write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Compute HF energy
ET = trace_matrix(nBas,matmul(P,T))
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EK = trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF)
end subroutine SPHF

71
src/QuAcK/SPMP2.f90
View File

@ -1,71 +0,0 @@
subroutine SPMP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform third-order Moller-Plesset calculation
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: ENuc,EHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps,E2a,E2b
! Output variables
double precision,intent(out) :: EcMP2(3)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset second-order calculation |'
write(*,*)'************************************************'
write(*,*)
! Compute MP2 energy
E2a = 0d0
E2b = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(i) + e(j) - e(a) - e(b)
! Secon-order ring diagram
E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps
! Second-order exchange
E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps
enddo
enddo
enddo
enddo
EcMP2(2) = E2a
EcMP2(3) = -E2b
EcMP2(1) = EcMP2(2) + EcMP2(3)
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part ',EcMP2(3)
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy',ENuc + EHF + EcMP2(1)
write(*,'(A32)') '-----------------------'
write(*,*)
end subroutine SPMP2

77
src/QuAcK/SPTDHF.f90
View File

@ -1,77 +0,0 @@
subroutine SPTDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI,e)
! Perform random phase approximation calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
logical :: dRPA,TDA,BSE
integer :: ispin
double precision,allocatable :: Omega(:,:),XpY(:,:,:)
double precision :: rho
double precision :: EcRPA
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Time-dependent Hartree-Fock calculation |'
write(*,*)'************************************************'
write(*,*)
! Switch on exchange for TDHF
dRPA = .false.
! Switch off Tamm-Dancoff approximation for TDHF
TDA = .false.
! Switch off Bethe-Salpeter equation for TDHF
BSE = .false.
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
end subroutine SPTDHF

81
src/QuAcK/SP_linear_response.f90
View File

@ -1,81 +0,0 @@
subroutine SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRPA,Omega,XpY)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS)
! Local variables
double precision :: trace_matrix
double precision,allocatable :: A(:,:),B(:,:),ApB(:,:),AmB(:,:),AmBSq(:,:),Z(:,:)
! Output variables
double precision,intent(out) :: EcRPA
double precision,intent(out) :: Omega(nS),XpY(nS,nS)
! Memory allocation
allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),Z(nS,nS))
! Build A and B matrices
call SP_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A)
! Tamm-Dancoff approximation
B = 0d0
if(.not. TDA) then
call SP_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B)
endif
! Build A + B and A - B matrices
AmB = A - B
ApB = A + B
! print*,'A+B'
! call matout(nS,nS,ApB)
! print*,'A-B'
! call matout(nS,nS,AmB)
! Diagonalize TD-HF matrix
call diagonalize_matrix(nS,AmB,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
call ADAt(nS,AmB,sqrt(Omega),AmBSq)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
call diagonalize_matrix(nS,Z,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
Omega = sqrt(Omega)
XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(Omega),XpY)
! print*,'RPA excitations'
! call matout(nS,1,Omega)
! Compute the RPA correlation energy
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
end subroutine SP_linear_response

56
src/QuAcK/SP_linear_response_A_matrix.f90
View File

@ -1,56 +0,0 @@
subroutine SP_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: Kronecker_delta
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: A_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ 0.5d0*(1d0 + delta_spin)*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*ERI(i,b,j,a)
enddo
enddo
enddo
enddo
end subroutine SP_linear_response_A_matrix

54
src/QuAcK/SP_linear_response_B_matrix.f90
View File

@ -1,54 +0,0 @@
subroutine SP_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: B_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
B_lr(ia,jb) = 0.5d0*(1d0 + delta_spin)*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*ERI(i,j,b,a)
enddo
enddo
enddo
enddo
end subroutine SP_linear_response_B_matrix