added sorted eigenvalues complex diag

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:")

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

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
+  complex*16                    :: ET
-  double precision              :: EV
+  complex*16                    :: EV
-  double precision              :: EJ
+  complex*16                    :: EJ
-  double precision              :: EK
+  complex*16                    :: EK
-  double precision              :: dipole(ncart)
+  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,Fp,eHF)
-    call complex_diagonalize_matrix(nBas,cp,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(*,*)'--------------------------------------------------------------------------------------------------'

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        :: cp(:,:)
-  complex*16,allocatable  :: e(:)
+  complex*16,allocatable        :: e(:)
 
 ! Output variables
 
-  complex*16,intent(out)  :: c(nBas,nOrb)
+  complex*16,intent(inout)  :: c(nBas,nOrb)
 
 ! Memory allocation
 

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
 

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)

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

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 

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

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