From a8dba8205c309a323d01bff21402fa9fc4d81196 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 19 Mar 2019 10:15:55 +0100 Subject: [PATCH] removing SP files --- src/QuAcK/Makefile | 36 +++++ src/QuAcK/SPHF.f90 | 170 ---------------------- src/QuAcK/SPMP2.f90 | 71 --------- src/QuAcK/SPTDHF.f90 | 77 ---------- src/QuAcK/SP_linear_response.f90 | 81 ----------- src/QuAcK/SP_linear_response_A_matrix.f90 | 56 ------- src/QuAcK/SP_linear_response_B_matrix.f90 | 54 ------- 7 files changed, 36 insertions(+), 509 deletions(-) create mode 100644 src/QuAcK/Makefile delete mode 100644 src/QuAcK/SPHF.f90 delete mode 100644 src/QuAcK/SPMP2.f90 delete mode 100644 src/QuAcK/SPTDHF.f90 delete mode 100644 src/QuAcK/SP_linear_response.f90 delete mode 100644 src/QuAcK/SP_linear_response_A_matrix.f90 delete mode 100644 src/QuAcK/SP_linear_response_B_matrix.f90 diff --git a/src/QuAcK/Makefile b/src/QuAcK/Makefile new file mode 100644 index 0000000..2c54de3 --- /dev/null +++ b/src/QuAcK/Makefile @@ -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 diff --git a/src/QuAcK/SPHF.f90 b/src/QuAcK/SPHF.f90 deleted file mode 100644 index e6cd623..0000000 --- a/src/QuAcK/SPHF.f90 +++ /dev/null @@ -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 diff --git a/src/QuAcK/SPMP2.f90 b/src/QuAcK/SPMP2.f90 deleted file mode 100644 index d91d803..0000000 --- a/src/QuAcK/SPMP2.f90 +++ /dev/null @@ -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 diff --git a/src/QuAcK/SPTDHF.f90 b/src/QuAcK/SPTDHF.f90 deleted file mode 100644 index 6409607..0000000 --- a/src/QuAcK/SPTDHF.f90 +++ /dev/null @@ -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 diff --git a/src/QuAcK/SP_linear_response.f90 b/src/QuAcK/SP_linear_response.f90 deleted file mode 100644 index e087397..0000000 --- a/src/QuAcK/SP_linear_response.f90 +++ /dev/null @@ -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 diff --git a/src/QuAcK/SP_linear_response_A_matrix.f90 b/src/QuAcK/SP_linear_response_A_matrix.f90 deleted file mode 100644 index d95ebf4..0000000 --- a/src/QuAcK/SP_linear_response_A_matrix.f90 +++ /dev/null @@ -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 diff --git a/src/QuAcK/SP_linear_response_B_matrix.f90 b/src/QuAcK/SP_linear_response_B_matrix.f90 deleted file mode 100644 index 6e60338..0000000 --- a/src/QuAcK/SP_linear_response_B_matrix.f90 +++ /dev/null @@ -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