diff --git a/crhf_test/generate_random_cap_dat.py b/crhf_test/generate_random_cap_dat.py index 64ae009..a38ce21 100644 --- a/crhf_test/generate_random_cap_dat.py +++ b/crhf_test/generate_random_cap_dat.py @@ -34,7 +34,7 @@ if __name__ == "__main__": with open("../int/nBas.dat", 'r') as f: size = int(f.readline().strip()) print("nBas: ", size) - width = 5 + width = 0 generate_cap_file("../int/CAP.dat", size, width) W = read_and_construct_matrix("../int/CAP.dat", size) print("W matrix:") diff --git a/crhf_test/run_test.sh b/crhf_test/run_test.sh index c44fa79..9c5b603 100755 --- a/crhf_test/run_test.sh +++ b/crhf_test/run_test.sh @@ -3,4 +3,4 @@ cp ./methods.test ../input/methods cp ./options.test ../input/options cd .. -python3 PyDuck.py -x N2 -b sto-3g -c -1 -m 2 +python3 PyDuck.py -x N2 -b sto-3g -c 0 -m 1 diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90 index 0e232cd..bc604a4 100644 --- a/src/HF/cRHF.f90 +++ b/src/HF/cRHF.f90 @@ -34,11 +34,11 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r integer :: nSCF integer :: nBasSq integer :: n_diis - double precision :: ET - double precision :: EV - double precision :: EJ - double precision :: EK - double precision :: dipole(ncart) + complex*16 :: ET + complex*16 :: EV + complex*16 :: EJ + complex*16 :: EK + complex*16 :: dipole(ncart) double precision :: Conv double precision :: rcond @@ -91,10 +91,10 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r Hc(:,:) = cmplx(T+V,W,kind=8) ! Guess coefficients and density matrix + call complex_mo_guess(nBas,nBas,guess_type,S,Hc,X,c) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - ! Initialization n_diis = 0 @@ -135,18 +135,21 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Kinetic energy - ET = trace_matrix(nBas,matmul(P,T)) + ET = cmplx(trace_matrix(nBas,real(matmul(P,T))),trace_matrix(nBas,aimag(matmul(P,T))),kind=8) ! Potential energy - EV = trace_matrix(nBas,matmul(P,V)) + EV = cmplx(trace_matrix(nBas,real(matmul(P,V))),trace_matrix(nBas,aimag(matmul(P,V))),kind=8) + ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + EJ = 0.5d0*cmplx(trace_matrix(nBas,real(matmul(P,J))),trace_matrix(nBas,aimag(matmul(P,J))),kind=8) + ! Exchange energy - EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + EK = 0.25d0*cmplx(trace_matrix(nBas,real(matmul(P,K))),trace_matrix(nBas,aimag(matmul(P,K))),kind=8) + ! Total energy @@ -160,16 +163,15 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! call complex_DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) ! ! end if -! + ! Level shift - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) call complex_level_shifting(level_shift,nBas,nBas,nO,S,c,F) ! Diagonalize Fock matrix Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - write(*,*) nBas - call complex_diagonalize_matrix(nBas,cp,eHF) + call complex_diagonalize_matrix(nBas,Fp,eHF) c = matmul(X,cp) ! Density matrix @@ -178,7 +180,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F32.10,1X,A1,1X,F32.10,A1,1X,A1,1X,F32.10,1X,A1,1X,F32.10,1X,A1,1X,E10.2,1X,A1,1X)') & - '|',nSCF,'|',real(ERHF),'+',aimag(ERHF),'i','|',EJ,'|',EK,'|',Conv,'|' + '|',nSCF,'|',real(ERHF),'+',aimag(ERHF),'i','|',real(EJ),'|',real(EK),'|',Conv,'|' write(*,*) real(ERHF),'+',aimag(ERHF),'i' end do write(*,*)'--------------------------------------------------------------------------------------------------' diff --git a/src/HF/complex_core_guess.f90 b/src/HF/complex_core_guess.f90 index 80b3f04..622cb36 100644 --- a/src/HF/complex_core_guess.f90 +++ b/src/HF/complex_core_guess.f90 @@ -12,12 +12,12 @@ subroutine complex_core_guess(nBas, nOrb, Hc, X, c) ! Local variables - complex*16,allocatable :: cp(:,:) - complex*16,allocatable :: e(:) + complex*16,allocatable :: cp(:,:) + complex*16,allocatable :: e(:) ! Output variables - complex*16,intent(out) :: c(nBas,nOrb) + complex*16,intent(inout) :: c(nBas,nOrb) ! Memory allocation diff --git a/src/HF/complex_mo_guess.f90 b/src/HF/complex_mo_guess.f90 index 9203f82..34e5157 100644 --- a/src/HF/complex_mo_guess.f90 +++ b/src/HF/complex_mo_guess.f90 @@ -28,9 +28,7 @@ subroutine complex_mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) write(*,*) 'Core guess...' call complex_core_guess(nBas, nOrb, Hc, X, c) - else - print*,'Wrong guess option' stop diff --git a/src/utils/complex_level_shifting.f90 b/src/utils/complex_level_shifting.f90 index 35b96ce..0bb478f 100644 --- a/src/utils/complex_level_shifting.f90 +++ b/src/utils/complex_level_shifting.f90 @@ -1,4 +1,4 @@ -subroutine level_shifting(level_shift, nBas, nOrb, nO, S, c, F) +subroutine complex_level_shifting(level_shift, nBas, nOrb, nO, S, c, F) ! Perform level-shifting on the Fock matrix @@ -21,7 +21,7 @@ subroutine level_shifting(level_shift, nBas, nOrb, nO, S, c, F) ! Output variables - double precision,intent(inout):: F(nBas,nBas) + complex*16,intent(inout) :: F(nBas,nBas) allocate(F_MO(nOrb,nOrb), Sc(nBas,nOrb)) complex_level_shift = cmplx(level_shift, 0.0,kind=8) diff --git a/src/utils/complex_matout.f90 b/src/utils/complex_matout.f90 new file mode 100644 index 0000000..a7a3dde --- /dev/null +++ b/src/utils/complex_matout.f90 @@ -0,0 +1,15 @@ +subroutine complex_matout(n,m,A) + + implicit none + +! Input variables + + integer,intent(in) :: n,m + complex*16,intent(in) :: A(n,m) + + write( *,*) 'Real part:' + call matout(n,m,real(A)) + write (*,*) 'Imaginary part:' + call matout(n,m,aimag(A)) + +end subroutine diff --git a/src/utils/complex_sort_eigenvalues.f90 b/src/utils/complex_sort_eigenvalues.f90 new file mode 100644 index 0000000..ae7f811 --- /dev/null +++ b/src/utils/complex_sort_eigenvalues.f90 @@ -0,0 +1,33 @@ +subroutine complex_sort_eigenvalues(N, eigvals, eigvecs) + implicit none + integer, intent(in) :: N + complex*16, intent(inout) :: eigvals(N) + complex*16, intent(inout) :: eigvecs(N, N) + + integer :: i, j + complex*16 :: temp_eigval + complex*16 :: temp_vec(N) + logical :: swapped + + do i = 1, N-1 + swapped = .FALSE. + do j = 1, N-i + if (REAL(eigvals(j)) > REAL(eigvals(j+1))) then + ! Swap eigenvalues + temp_eigval = eigvals(j) + eigvals(j) = eigvals(j+1) + eigvals(j+1) = temp_eigval + + ! Swap corresponding eigenvectors + temp_vec = eigvecs(:, j) + eigvecs(:, j) = eigvecs(:, j+1) + eigvecs(:, j+1) = temp_vec + + swapped = .TRUE. + end if + end do + ! If no swaps were made, the array is already sorted + if (.not. swapped) exit + end do + +end subroutine diff --git a/src/utils/complex_vecout.f90 b/src/utils/complex_vecout.f90 new file mode 100644 index 0000000..d910ebf --- /dev/null +++ b/src/utils/complex_vecout.f90 @@ -0,0 +1,17 @@ +subroutine complex_vecout(n,v) + + implicit none + +! Input variables + integer,intent(in) :: n + complex*16,intent(in) :: v(n) +! Local variables + double precision,allocatable :: v2(:,:) + + allocate(v2(n,2)) + write(*,*) 'First column real part, second imaginary part' + v2(:,1) = real(v) + v2(:,2) = aimag(v) + call matout(n,2,v2) + deallocate(v2) +end subroutine diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 5077069..17b5c18 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -149,23 +149,22 @@ subroutine complex_diagonalize_matrix(N,A,e) complex*16 :: VR(N,N) ! Memory allocation - allocate(work(1)) - lwork = -1 call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info) - lwork = int(work(1)) + lwork = max(1,int(real(work(1)))) + deallocate(work) - allocate(work(lwork)) call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info) + call complex_sort_eigenvalues(N,e,VR) deallocate(work) A = VR if(info /= 0) then - print*,'Problem in diagonalize_matrix (zgees)!!' + print*,'Problem in diagonalize_matrix (zgeev)!!' end if end subroutine