PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:24:58 +02:00
Code Releases Activity

First debugged version of cRHF

Browse Source
This commit is contained in:
Loris Burth 2025-03-10 14:46:56 +01:00
parent 06ef322a39
commit 20fe1bc59a
14 changed files with 382 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

209
crhf_test/cRHF.py Executable file
View File

@ -0,0 +1,209 @@
#!/usr/bin/env python3
import numpy as np
import pandas as pd
import os
def print_matrix(matrix):
df = pd.DataFrame(matrix)
print(df)
def read_nO(working_dir):
file_path = os.path.join(working_dir, "input/molecule")
with open(file_path, "r") as f:
next(f) # Skip the first line
line = f.readline().split()
nO = max(int(line[1]), int(line[2]))
return nO
def read_ENuc(file_path):
# Path to the ENuc.dat file
with open(file_path, 'r') as f:
# Read the nuclear repulsion energy from the first line
ENuc = float(f.readline().strip())
return ENuc
def read_matrix(filename):
# Read the data and determine matrix size
entries = []
max_index = 0
with open(filename, 'r') as f:
for line in f:
i, j, value = line.split()
i, j = int(i) - 1, int(j) - 1 # Convert to zero-based index
entries.append((i, j, float(value)))
# Find max index to determine size
max_index = max(max_index, i, j)
# Initialize square matrix with zeros
matrix = np.zeros((max_index + 1, max_index + 1))
# Fill the matrix
for i, j, value in entries:
matrix[i, j] = value
if i != j: # Assuming the matrix is symmetric, fill the transpose element
matrix[j, i] = value
return matrix
def read_CAP_integrals(filename, size):
"""
Reads the file and constructs the symmetric matrix W.
"""
W = np.zeros((size, size))
with open(filename, 'r') as f:
for line in f:
mu, nu, wx, wy, wz = line.split()
mu, nu = int(mu) - 1, int(nu) - 1 # Convert to zero-based index
value = float(wx) + float(wy) + float(wz)
W[mu, nu] = value
W[nu, mu] = value # Enforce symmetry
return W
def read_2e_integrals(file_path, nBas):
# Read the binary file and reshape the data into a 4D array
try:
G = np.fromfile(file_path, dtype=np.float64).reshape(
(nBas, nBas, nBas, nBas))
except FileNotFoundError:
print(f"Error opening file: {file_path}")
raise
return G
def get_X(S):
"""
Computes matrix X for orthogonalization. Attention O has to be hermitian.
"""
vals, U = np.linalg.eigh(S)
# Sort the eigenvalues and eigenvectors
vals = 1/np.sqrt(vals)
return U@np.diag(vals)
def sort_eigenpairs(eigenvalues, eigenvectors):
# Get the sorting order based on the real part of the eigenvalues
order = np.argsort(eigenvalues.real)
# Sort eigenvalues and eigenvectors
sorted_eigenvalues = eigenvalues[order]
sorted_eigenvectors = eigenvectors[:, order]
return sorted_eigenvalues, sorted_eigenvectors
def diagonalize(M):
# Diagonalize the matrix
vals, vecs = np.linalg.eig(M)
# Sort the eigenvalues and eigenvectors
vals, vecs = sort_eigenpairs(vals, vecs)
# Orthonormalize them wrt cTc inner product
vecs = gram_schmidt(vecs)
return vals, vecs
def Hartree_matrix_AO_basis(P, ERI):
# Initialize Hartree matrix with zeros (complex type)
J = np.zeros((nBas, nBas), dtype=np.complex128)
# Compute Hartree matrix
for si in range(nBas):
for nu in range(nBas):
for la in range(nBas):
for mu in range(nBas):
J[mu, nu] += P[la, si] * ERI[mu, la, nu, si]
return J
def exchange_matrix_AO_basis(P, ERI):
# Initialize exchange matrix with zeros
K = np.zeros((nBas, nBas), dtype=np.complex128)
# Compute exchange matrix
for nu in range(nBas):
for si in range(nBas):
for la in range(nBas):
for mu in range(nBas):
K[mu, nu] -= P[la, si] * ERI[mu, la, si, nu]
return K
def gram_schmidt(vectors):
"""
Orthonormalize a set of vectors with respect to the scalar product c^T c.
"""
orthonormal_basis = []
for v in vectors.T: # Iterate over column vectors
for u in orthonormal_basis:
v -= (u.T @ v) * u # Projection with respect to c^T c
norm = np.sqrt(v.T @ v) # Norm with respect to c^T c
if norm > 1e-10:
orthonormal_basis.append(v / norm)
else:
raise Exception("Norm of eigenvector < 1e-10")
return np.column_stack(orthonormal_basis)
if __name__ == "__main__":
# Constants
workdir = "../"
eta = 0.01
thresh = 0.00001
maxSCF = 256
nO = read_nO(workdir)
# Read integrals
T = read_matrix("../int/Kin.dat")
S = read_matrix("../int/Ov.dat")
V = read_matrix("../int/Nuc.dat")
ENuc = read_ENuc("../int/ENuc.dat")
nBas = np.shape(T)[0]
W = read_CAP_integrals("../int/CAP.dat", nBas)
ERI = read_2e_integrals("../int/ERI.bin", nBas)
X = get_X(S)
W = -eta*W
Hc = T + V + W*1j
# core guess
_, c = diagonalize(X.T @ Hc @ X)
c = X @ c
P = 2*c[:, :nO]@c[:, :nO].T
print('-' * 98)
print(
f"| {'#':<1} | {'E(RHF)':<36} | {'EJ(RHF)':<16} | {'EK(RHF)':<16} | {'Conv':<10} |")
print('-' * 98)
nSCF = 0
Conv = 1
while(Conv > thresh and nSCF < maxSCF):
nSCF += 1
J = Hartree_matrix_AO_basis(P, ERI)
K = exchange_matrix_AO_basis(P, ERI)
F = Hc + J + 0.5*K
err = F@P@S - S@P@F
if nSCF > 1:
Conv = np.max(np.abs(err))
ET = np.trace(P@T)
EV = np.trace(P@V)
EJ = 0.5*np.trace(P@J)
EK = 0.25*np.trace(P@K)
ERHF = ET + EV + EJ + EK
Fp = X.T @ F @ X
eHF, c = diagonalize(Fp)
c = X @ c
P = 2*c[:, :nO]@c[:, :nO].T
print(
f"| {nSCF:3d} | {ERHF.real+ENuc:5.6f} + {ERHF.imag:5.6f}i | {EJ:5.6f} | {EK:5.6f} | {Conv:5.6f} |")
print('-' * 98)
print()
print("RHF orbitals")
print_matrix(eHF)

1
crhf_test/matrix_test
View File

File diff suppressed because one or more lines are too long

36
crhf_test/test_diagonalization.py
View File

@ -1,36 +0,0 @@
import numpy as np
import re
def read_complex_matrix(filename, N):
with open(filename, "r") as file:
lines = file.readlines()
# Flattened list of complex numbers
complex_numbers = []
# Regular expression to match complex numbers in the format (a,b)
complex_pattern = re.compile(
r"\(\s*([-+]?\d*\.?\d+(?:[Ee][-+]?\d+)?),\s*([-+]?\d*\.?\d+(?:[Ee][-+]?\d+)?)\)")
for line in lines:
matches = complex_pattern.findall(line)
for real, imag in matches:
complex_numbers.append(complex(float(real), float(imag)))
# Reshape into NxN matrix
return np.array(complex_numbers, dtype=np.complex128).reshape(N, N)
# Read the matrix stored from the first core guess and compare values to the fortran code --> wokrs
N = 10 # Set matrix size
matrix = read_complex_matrix("matrix_test", N)
print("Matrix")
print(matrix)
print("Element 5,1")
print(matrix[4, 0])
eigenval, eigenvec = np.linalg.eig(matrix)
print("Eigenvalues")
print(eigenval)
print("Eigenvectors")
print(eigenvec)

3
src/AOtoMO/complex_exchange_matrix_AO_basis.f90
View File

@ -25,8 +25,7 @@ subroutine complex_exchange_matrix_AO_basis(nBas,P,ERI,K)
do si=1,nBas
do la=1,nBas
do mu=1,nBas
K(mu,nu) = K(mu,nu) - P(la,si)
!K(mu,nu) = K(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
K(mu,nu) = K(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
end do
end do
end do

7
src/HF/RHF.f90
View File

@ -89,10 +89,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Guess coefficients and density matrix
call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
write(*,*)'Initial guess c'
call matout(nBas,nBas,c)
write(*,*) "verify orthonormalization"
call matout(nBas,nBas,matmul(transpose(c),matmul(S,c)))
P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
@ -131,8 +127,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
write(*,*)'Fock matrix'
call matout(nBas,nBas,F)
! Check convergence
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)

50
src/HF/cRHF.f90
View File

@ -29,7 +29,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
! Local variables
integer :: nSCF
integer :: nBasSq
@ -82,33 +82,18 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis), &
Hc(nBas,nBas),W(nBas,nBas),Z(nBas,nBas))
! Read CAP integrals from file
call read_CAP_integrals(nBas,W)
W(:,:) = -eta*W(:,:)
! Define core Hamiltonian with CAP part
Hc(:,:) = cmplx(T+V,W,kind=8)
! Guess coefficients and density matrix
call complex_mo_guess(nBas,nBas,guess_type,S,Hc,X,c)
write(*,*) "Guess coefficients c"
call complex_matout(nBas,nBas,c)
write(*,*) "Check if c orthonormal"
call complex_matout(nBas,nBas,matmul(transpose(c),matmul(S,c)))
call complex_orthogonalize_matrix(nBas,matmul(transpose(c),matmul(S,c)),Z)
c = matmul(c,Z)
write(*,*) "Check if c tilde orthonormal"
call complex_matout(nBas,nBas,matmul(transpose(c),matmul(S,c)))
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
! Initialization
n_diis = 0
F_diis(:,:) = cmplx(0d0,0d0,kind=8)
err_diis(:,:) = cmplx(0d0,0d0,kind=8)
@ -116,7 +101,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
@ -133,13 +117,12 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
nSCF = nSCF + 1
! Build Fock matrix
call complex_Hartree_matrix_AO_basis(nBas,P,ERI,J)
call complex_exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
write(*,*) "Fock matrix"
call complex_matout(nBas,nBas,F)
! Check convergence
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
@ -167,34 +150,31 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
ERHF = ET + EV + EJ + EK
! ! DIIS extrapolation (fix later) !
!
! DIIS extrapolation (fix later) !
! if(max_diis > 1) then
!
! n_diis = min(n_diis+1,max_diis)
! 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 complex_level_shifting(level_shift,nBas,nBas,nO,S,c,F)
!
! ! Level shift
! 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))
Fp = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
write(*,*)'t1'
call complex_diagonalize_matrix(nBas,Fp,eHF)
call complex_diagonalize_matrix(nBas,cp,eHF)
c = matmul(X,cp)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! 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','|',real(EJ),'|',real(EK),'|',Conv,'|'
write(*,*) real(ERHF),'+',aimag(ERHF),'i'
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,A1,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',real(ERHF + ENuc),'+',aimag(ERHF),'i','|',real(EJ),'|',real(EK),'|',Conv,'|'
end do
write(*,*)'--------------------------------------------------------------------------------------------------'
!------------------------------------------------------------------------
@ -217,7 +197,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call complex_dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Testing zone

53
src/HF/complex_dipole_moment.f90 Normal file
View File

@ -0,0 +1,53 @@
subroutine complex_dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
! Compute density matrix based on the occupation numbers
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nNuc
complex*16,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ixyz
integer :: iNuc
integer :: mu,nu
! Output variables
complex*16,intent(out) :: dipole(ncart)
! Initialization
dipole(:) = 0d0
! Loop over cartesian components
do ixyz=1,ncart
! Nuclear part
do iNuc=1,nNuc
dipole(ixyz) = dipole(ixyz) + ZNuc(iNuc)*rNuc(iNuc,ixyz)
end do
! Electronic part
do nu=1,nBas
do mu=1,nBas
dipole(ixyz) = dipole(ixyz) - P(mu,nu)*dipole_int(mu,nu,ixyz)
end do
end do
end do
end subroutine

4
src/HF/core_guess.f90
View File

@ -29,6 +29,10 @@ subroutine core_guess(nBas, nOrb, Hc, X, c)
cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:)))
call diagonalize_matrix(nOrb, cp, e)
write(*,*) 'cp'
call matout(nOrb,nOrb,cp)
write(*,*) 'Eigenvalues e'
call vecout(nOrb,e)
c(:,:) = matmul(X(:,:), cp(:,:))
deallocate(cp, e)

60
src/HF/print_cRHF.f90
View File

@ -12,14 +12,14 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
complex*16,intent(in) :: eHF(nOrb)
complex*16,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
double precision,intent(in) :: EJ
double precision,intent(in) :: EK
double precision,intent(in) :: ERHF
complex*16,intent(in) :: ET
complex*16,intent(in) :: EJ
complex*16,intent(in) :: EK
complex*16,intent(in) :: EV
complex*16,intent(in) :: ERHF
double precision,intent(in) :: dipole(ncart)
! Local variables
@ -27,7 +27,7 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
integer :: ixyz
integer :: HOMO
integer :: LUMO
double precision :: Gap
complex*16 :: Gap
double precision :: S,S2
logical :: dump_orb = .false.
@ -44,32 +44,32 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A33)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au'
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au'
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au'
write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au'
write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',ERHF,' au'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV),'+',aimag(ET+EV),' au'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Kinetic energy = ',real(ET),'+',aimag(ET),' au'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Potential energy = ',real(EV),'+',aimag(Ev),'i',' au'
write(*,'(A50)') '------------------------------------------------------'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Two-electron energy = ',real(EJ + EK),'+',aimag(EJ+EK),'i',' au'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Hartree energy = ',real(EJ),'+',aimag(EJ),'i',' au'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Exchange energy = ',real(EK),'+',aimag(EK),'i',' au'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' Electronic energy = ',real(ERHF),'+',aimag(ERHF),'i',' au'
write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
write(*,'(A33,1X,F16.10,A3)') ' cRHF energy = ',ERHF + ENuc,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6,A3)') ' HF HOMO energy = ',eHF(HOMO)*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HF LUMO energy = ',eHF(LUMO)*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HF HOMO-LUMO gap = ',Gap*HaToeV,' eV'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A1,F16.10,A1,A3)') ' cRHF energy = ',real(ERHF + ENuc),'+',aimag(ERHF+ENuc),'i',' au'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A33,1X,F16.6,A3)') ' HF HOMO energy = ',real(eHF(HOMO))*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HF LUMO energy = ',real(eHF(LUMO))*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HF HOMO-LUMO gap = ',real(Gap)*HaToeV,' eV'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A33,1X,F16.6)') ' <Sz> = ',S
write(*,'(A33,1X,F16.6)') ' <S^2> = ',S2
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') '------------------------------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '---------------------------------------'
write(*,'(10X,4F10.4)') (real(dipole(ixyz))*auToD,ixyz=1,ncart),norm2(real(dipole))*auToD
write(*,'(A50)') '---------------------------------------------'
write(*,*)
! Print results
@ -78,13 +78,13 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' cRHF orbital coefficients '
write(*,'(A50)') '---------------------------------------'
call matout(nBas, nOrb, cHF)
call complex_matout(nBas, nOrb, cHF)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' cRHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nOrb, eHF)
call complex_vecout(nOrb, eHF)
write(*,*)
end subroutine

14
src/QuAcK/RQuAcK.f90
View File

@ -98,10 +98,14 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:)
complex*16,allocatable :: complex_eHF(:)
double precision,allocatable :: cHF(:,:)
complex*16,allocatable :: complex_cHF(:,:)
double precision,allocatable :: PHF(:,:)
complex*16,allocatable :: complex_PHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision :: ERHF
complex*16 :: complex_ERHF
double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
@ -117,15 +121,16 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
!-------------------!
! Memory allocation !
!-------------------!
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas))
allocate(complex_eHF(nOrb))
allocate(complex_cHF(nBas,nOrb))
allocate(complex_PHF(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
allocate(FHF(nBas,nBas))
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
@ -165,10 +170,9 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
end if
if(docRHF) then
call wall_time(start_HF)
call cRHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
nBas,nO,S,T,V,ERI_AO,dipole_int_AO,X,ERHF,complex_eHF,complex_cHF,complex_PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF

55
src/utils/complex_gram_schmidt.f90 Normal file
View File

@ -0,0 +1,55 @@
subroutine complex_gram_schmidt(N, vectors)
! Input variables
implicit none
integer, intent(in) :: N
complex*16, intent(inout) :: vectors(N, N)
! Local variables
integer :: i, j
complex*16 :: Mtmp(N,N)
complex*16 :: proj
complex*16 :: norm
! Copy the input matrix to a temporary matrix
Mtmp(:, :) = vectors(:,:)
! Orthonormalize the vectors
do i = 1, N
! Orthonormalize with respect to the previous vectors in Mtmp
do j = 1, i-1
! Compute the dot product (scalar product) of vectors j and i
call complex_dot_product(N,Mtmp(:, j), Mtmp(:, i),proj)
! Subtract the projection onto the j-th vector
Mtmp(:, i) = Mtmp(:, i) - proj * Mtmp(:, j)
end do
! Normalize the vector
call complex_dot_product(N,Mtmp(:, i), Mtmp(:, i),proj)
norm = sqrt(proj)
if (abs(norm) > 1.0d-10) then
! Normalize the i-th vector and store it back in vectors
vectors(:, i) = Mtmp(:, i) / norm
else
print*, "Error: Norm of eigenvector < 1e-10 !!!"
stop
end if
end do
end subroutine
subroutine complex_dot_product(N,vec1,vec2,proj)
! Input
complex*16,intent(in) :: vec1(N),vec2(N)
!Output
complex*16, intent(out) :: proj
! Local variables
integer :: i
proj = 0d0
do i=1,N
proj = proj +vec1(i)*vec2(i)
end do
end subroutine

14
src/utils/complex_orthogonalize_matrix.f90
View File

@ -26,23 +26,11 @@ subroutine complex_orthogonalize_matrix(nBas,S,X)
allocate(Uvec(nBas,nBas),Uval(nBas),D(nBas,nBas))
! write(*,*)
! write(*,*) ' Lowdin orthogonalization'
! write(*,*)
Uvec = S
call complex_diagonalize_matrix(nBas,Uvec,Uval)
call complex_matout(nBas,nBas,matmul(Uvec,transpose(Uvec)))
do i=1,nBas
if(abs(Uval(i)) < thresh) then
write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i)
endif
Uval(i) = 1d0/sqrt(Uval(i))
enddo
call diag(nBas,Uval, D)
X = matmul(Uvec,matmul(D,transpose(Uvec)))

1
src/utils/complex_sort_eigenvalues.f90
View File

@ -1,4 +1,5 @@
subroutine complex_sort_eigenvalues(N, eigvals, eigvecs)
! Sort eigenvalues and corresponding eigenvectors wrt the real part of the eigenvalue
implicit none
integer, intent(in) :: N
complex*16, intent(inout) :: eigvals(N)

4
src/utils/wrap_lapack.f90
View File

@ -159,7 +159,7 @@ subroutine complex_diagonalize_matrix(N,A,e)
call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info)
call complex_sort_eigenvalues(N,e,VR)
call complex_gram_schmidt(N,VR)
deallocate(work)
A = VR
@ -169,6 +169,8 @@ subroutine complex_diagonalize_matrix(N,A,e)
end subroutine
subroutine svd(N,A,U,D,Vt)
! Compute A = U.D.Vt